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

      Real*8, Allocatable :: xdat(:,:), ydat(:,:), xtb(:), sx(:)
      Real*8, Allocatable :: EOFs(:,:), PCs(:,:)
      Real*8, Allocatable   :: Eigenvalues(:), ECumulate(:)
      !---------------------------------------------------------------------------------

      Integer :: Yr1, Yr2, Nmon, Nyear ! Thoi ky
      Integer :: Nt, Nx   ! Kich thuoc thuc cua mang X
      Integer :: Mcount, Ncount, switch, Standard, PC_Stard, Nmode, K
      Character (Len=250) :: Fname, Fname1, MaskFname, pca_out, eof_out
      Real*8 :: Factor

      Namelist /Model_Info/ Fname, MaskFname, pca_out, eof_out, &
                            Yr1, Yr2, Nmon, &
                            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"/)
      !---------------------------------------------------------------------------------

      Open (1,file="Namelist_PCA_New.txt",Status="old")
      Read (1, Model_Info)
      Close (1)
      !---------------------------------------------------------------------------------
      ! Reading data:
      Call rd_ncdf_2d (trim(MaskFname),xlon,xlat,Nlon,Nlat,Mask, "land")
      Call Read_Drough_IDX (Fname)
      Nyear = Yr2-Yr1+1            ! So nam co SL
      Nt = Ntime; Nx = Nst
      !---------------------------------------------------------------------------------
      if (switch==1) then
         Mcount=Nt; Ncount=Nx
      else if (switch==0) then
         Mcount=Nx; Ncount=Nt
      endif
      Factor = Sqrt(Real(Nt))

      Allocate (xdat(Nt,Nx), ydat(Nt,Nx))
      Allocate (Xtb(Nx), Sx(Nx))   ! TB, DLC, ma tran Covar hoac Corre
      Allocate (Eigenvalues(Mcount))          ! Tri rieng va Vecto rieng (1D)
      Allocate (EOFs(Nt,Nmode), PCs(Nx,Nmode), ECumulate(Nmode))
      !---------------------------------------------------------------------------------
      ! Loop over all drought indices:
      !Fname="pca_dr_APH.nc"; Fname1="eof_timeseries_dr_APH.nc"
      Do iidx=1,NDr_idx
         Xdat(1:Nt,1:Nx) = Dr_IDX(1:Nt,1:Nx,iidx)
         ! Tinh cac EOFs, PCs
         Call PCA_2019 (Xdat,Ydat,Nt,Nx,switch, Standard, PC_Stard, &
                        Xtb,Sx,Eigenvalues,EOFs,PCs,Nmode) 

         Call Write_EOF_PCA (pca_out,eof_out,iidx)
         print*, "===========> Complete for ",IDX_nm(iidx)
      Enddo
      !Deallocate(Xbar,Std,Rx,EValues,V,Fx,Ax,D,EVtor,Lamda,EOF,PCs,ECumulate)
      print*,"Program is normal terminating ..."

      CONTAINS
      !================================================================================= 
      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,it
         Real :: tmp 
         Real, Allocatable :: PCA(:,:,:), xmode(:), EOF1(:,:), X1(:), Xtmp(:)

         Allocate (PCA(Nlon,Nlat,Nmode), xmode(Nmode), Xtmp(Nmode)) !, EOF1(Nmode,Nt))
         Allocate (EOF1(Nt,Nmode))
         !-------------------------------------------------------------------------         
         ! Write to PCs file:
         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
         !-------------------------------------------------------------------------
         Xtmp(1:Nmode) = Eigenvalues(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                 trim(IDX_nm(iidx))//"_eig", Nmode, &
                "Eigen values of "//trim(IDX_nm(iidx)), "-", "mode", Miss)

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

         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                trim(IDX_nm(iidx))//"_var",Nmode, &
                "Percent Contribution of "//trim(IDX_nm(iidx)), &
                "-", "mode", Miss )

         Xtmp(1:Nmode) = ECumulate(1:Nmode)
         Call wr_ncdf_1d_apd (F_pca, Xtmp, &
                trim(IDX_nm(iidx))//"_cum",Nmode, &
                "Cumuluted Eigenvalues of "//trim(IDX_nm(iidx)), &
                "-", "mode", Miss )
         !-------------------------------------------------------------------------
         ! Write to EOF file:
         if (Allocated(xtime)) Deallocate (xtime); Allocate (xtime(1:Nt))
         Do it=1, Nt; xtime(it)=it; Enddo     ! Months
         !Do i=1, Nmode; EOF1(i,1:Nt)=EOFs(1:Nt,i); Enddo
         Do i=1, Nmode; EOF1(1:Nt,i)=EOFs(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)
            Call wr_ncdf_2d (F_eof,xtime,xmode, Nt,Nmode, EOF1, &
                     trim(IDX_nm(iidx)), &
                     trim(var_ln), "-", "mon","mode", Miss)
         else
            !Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nmode,Nt, &
            !         trim(var_ln), "-", "mode", "mon", Miss)
            Call wr_ncdf_2d_apd (F_eof,EOF1, trim(IDX_nm(iidx)), Nt,Nmode, &
                     trim(var_ln), "-", "mon", "mode", Miss)
         endif
         print*,"Write_EOF_PCA: Done !!!"
      End Subroutine Write_EOF_PCA
      !---------------------------------------------------------------------------------
      Subroutine Read_Drough_IDX (Fname)
         Integer :: ist, ilon, ilat, it, Nt, imon, imon1, iyr
         Integer :: nlat1, nlat0
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         Real, Allocatable :: X0(:), X1(:,:,:), X2(:,:,:), X3(:,:,:)
         Real, Allocatable :: xlon(:), xlat(:), xtime(:)

         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")
         !Call rd_ncdf_2d (trim(MaskFname), xlon, xlat, nlon,nlat,Mask,"land")

         Ntime = Nt - Nmon
         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)

         X1(:,1:nlat,1:Ntime)=PDSI  (:,1:Nlat,Nmon+1:Nt)
         X2(:,1:nlat,1:Ntime)=SPEI12(:,1:Nlat,Nmon+1:Nt)
         X3(:,1:nlat,1:Ntime)=SPI12 (:,1:Nlat,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
         print*,"Nst= ",Nst
         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

      !---------------------------------------------------------------------------------

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

      END