Use NetCDF
      Implicit None


      Real, Parameter :: Miss=-99.0
      integer, parameter :: nmodes = 14
      integer :: nt, nx, iminnxnt, imaxnxnt, nx_real, ny_real
      double precision, allocatable :: data(:,:), eigenvalues(:), &
                              eigenvectors(:,:), &
                              princomp(:,:), variance(:), &
                              cumvariance(:)
      double precision :: value, time
      double precision, allocatable :: spacked(:), &
                             evals(:), evecs(:,:), &
                             rlawork(:), tentpcs(:,:), &
                             sqrootweights(:), lon(:), &
                             lat(:)
      integer, allocatable :: ilawork(:), &
                 ifail(:)
      integer :: i, j, k, idx

      integer :: icovcor
!------------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:)
      Real, Allocatable :: Xmean(:,:), SS(:,:)
      Integer :: Nlon, Nlat, Ntime
      Character (Len=100) :: Fname, Fname1
      Logical :: Status
!------------------------------------------------------------------------------------
      Integer, Parameter :: Nidx=14, NYr1=1980, NYr2=2017
      Character (Len=6), Parameter :: IDXName(Nidx) = &
           (/"ONI","MEI","REQSOI", "SOI", "RINDO",    &
             "Nino34", "BEST",  "PDO",    "PNA", "DMI", &
             "WHWP",   "TNI",   "WP",     "QBO"/)
      Real, allocatable :: ClimIDX(:,:), ClimMean(:,:) !, xlon(:), xlat(:)
      Real, allocatable :: ClimCentral(:,:), ClimStandard(:,:)
!------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------
      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:)
      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"/)
!------------------------------------------------------------------------------------
      Fname="data_test.txt" !"/home/tanpv/DATA/Clim_IDX.nc" !"/home/tanpv/Progs/Dr_SEA/Results/spi.nc"
      !Call Read_Drough_IDX (Fname)
      !Call Read_Clim_IDX (Fname)
      Call Read_Data_test (Fname)
      print*,"Reading data ... done !"
      !------------------------------------------------------------------------------
      ! Identify the sizes of all arrays:
      Nt = Ntime; nx_real=Nlon; ny_real=Nlat
      Nx = nx_real*ny_real
      iminnxnt=min(nx,nt); imaxnxnt=max(nx,nt)

      ! For Dr_IDX:
      !Nx = Nst ! nx_real*ny_real
      !iminnxnt=min(nx,nt); imaxnxnt=max(nx,nt)


      ! Allocate:
      Allocate (data(nx,nt),eigenvalues(nmodes),eigenvectors(nx,nmodes), &
               princomp(nt,nmodes), variance(nmodes),cumvariance(nmodes))
      Allocate (spacked(iminnxnt*(iminnxnt+1)/2),evals(iminnxnt), &
               evecs(iminnxnt,nmodes),rlawork(8*iminnxnt), &
               tentpcs(imaxnxnt,nmodes), &
               sqrootweights(nx), lon(nx_real), lat(ny_real), &
               ilawork(5*iminnxnt), ifail(iminnxnt))
      !Allocate (SS(nx,nx))
      !------------------------------------------------------------------------------
      ! Loop over all drought indices:
      !Do iidx=1, NDr_idx
         !Data(1:Nx,1:Nt) = Dr_IDX(1:Nst,1:Ntime,iidx)
         Do i=1, Nx
            Data(i,1:Nt) = ClimCentral(1:Nt,i)
         Enddo

         do i=1, nx
            sqrootweights(i) = 1.0   ! Set area weights to 1.0
         enddo
         icovcor = 0       ! use covariance matrix
         !icovcor = 1       ! use correlation matrix

         Call deof( data, nx, nt, nmodes, icovcor,eigenvalues, &
                 eigenvectors, &
                 princomp, variance, cumvariance, &
                 iminnxnt, imaxnxnt, spacked, evals, evecs, rlawork, &
                 ilawork, ifail, tentpcs, sqrootweights )
         print*,"EOF calculating ... Done !!!"
      !------------------------------------------------------------------------------
          ! Write results to files:
         ! Fname="pca_tx.nc"; Fname1="eof_tx.nc"
          Fname="pca_xt.nc"; Fname1="eof_xt.nc"
          Call Write_EOF_PCA (Fname,Fname1, 1) !,iidx)
      !Enddo


      CONTAINS
!----------------------------------------------------------------------------------!
      Subroutine Read_Data_test (Fname)
      Character (Len=*) :: Fname
      Real, Allocatable :: X(:,:), tmp(:)
      Integer :: iidx, iyr, imon, year, IO, it, it0, it1

      Open (1, file=Fname, status="Old")
      Read (1,*) Ntime, Nlon
      Allocate (X(Ntime,Nlon)); Nlat=1
      Do it=1, Ntime
         Read (1,*) (X(it,iidx),iidx=1,Nlon)
      Enddo 
      if (.not.Allocated(ClimIDX)) Allocate (ClimIDX(Ntime,Nlon))
      ClimIDX(:,:)=X(:,:)
      Allocate(ClimCentral(Ntime,Nlon),ClimStandard(Ntime,Nlon))
      ClimCentral(:,:)=ClimIDX(:,:)
      Call Centralize (ClimCentral, Ntime, Nlon, ClimStandard, Miss)
      !-------------------------------------------------------------
      ! For writing to files:
      Allocate (xlon(Nlon), xlat(Nlat), xtime(Ntime), ClimMean(Nlon,Nlat))
      if (.not.Allocated(tmp)) Allocate (tmp(Nlon))
      !Allocate (lon(Nlon), lat(Nlat))
      Do iidx=1, Nlon
         xlon(iidx)=iidx; ClimMean(iidx,1)=0
         tmp(iidx) = Ave(ClimIDX(1:Ntime,iidx),Ntime,Miss)  ! Trung binh chuoi thoi gian
         Do it=1, Ntime
            ClimMean(iidx,1)=ClimMean(iidx,1)+ClimIDX(it,iidx)
         Enddo
         ClimMean(iidx,1)=ClimMean(iidx,1)/Real(Ntime)
         Do it=1, Ntime
         !   ClimIDX(it,iidx)=ClimIDX(it,iidx)-ClimMean(iidx,1)
            ClimIDX(it,iidx)=ClimIDX(it,iidx)-tmp(iidx)         ! Anomaly
         Enddo
      Enddo
      xlat(:)=1
      Lon(:)=xlon(:); Lat(:)=xlat(:)
      it1=0
      Do iyr=1, Ntime
         xtime(iyr)=iyr
      Enddo

      print*, Ntime, Nlon, Nlat
      End Subroutine Read_Data_test
!----------------------------------------------------------------------------------!
      Subroutine Read_Clim_IDX (Fname)
      ! Read data of Climate Indices from NetCDF and store in 3D array of ClimIDX(1:Nyear,1:Nmon,1:Ntr)
      ! where Ntr=14 is number of Climate indices
      Integer :: iidx, iyr, imon, year, IO, it, it0, it1
      Character (Len=*) :: Fname
      Real, Allocatable :: X(:,:),Xyrs(:),XMon(:), tmp(:)
      Integer :: NYrs, NMons

      if (.not.Allocated(tmp)) Allocate (tmp(Nidx))
      Do iidx=1, Nidx
         if (Allocated(X)) Deallocate (X)
         if (Allocated(Xyrs)) Deallocate (Xyrs)
         if (Allocated(XMon)) Deallocate (XMon)

         Call rd_ncdf_2d(trim(Fname),XYrs,XMon,NYrs,NMons, X, trim(IDXName(iidx)))
         if (.not.Allocated(ClimIDX)) Allocate (ClimIDX((NYr2-NYr1+1)*NMons,Nidx))
         it1=0
         Do iyr=NYr1, NYr2
            Do it=1, NYrs
               if (Int(XYrs(it))==iyr) then
                  Do imon=1, NMons
                     it1=it1+1; ClimIDX(it1,iidx)=X(it,imon)
                  Enddo
               endif
            Enddo
         Enddo
      Enddo
      Ntime = (NYr2-NYr1+1)*NMons; Nlon=Nidx; Nlat=1
      ! Qui tam va Chuan hoa: ClimCentral(:,:), ClimStandard(:,:)
      !-------------------------------------------------------------
      Allocate(ClimCentral(Ntime,Nidx),ClimStandard(Ntime,Nidx))
      ClimCentral(:,:)=ClimIDX(:,:)
      Call Centralize (ClimCentral, Ntime, Nidx, ClimStandard, Miss)
      !-------------------------------------------------------------
      ! For writing to files:
      Allocate (xlon(Nlon), xlat(Nlat), xtime(Ntime), ClimMean(Nlon,Nlat))
      !Allocate (lon(Nlon), lat(Nlat))
      Do iidx=1, Nidx
         xlon(iidx)=iidx; ClimMean(iidx,1)=0
         tmp(iidx) = Ave(ClimIDX(1:Ntime,iidx),Ntime,Miss)  ! Trung binh chuoi thoi gian
         Do it=1, Ntime
            ClimMean(iidx,1)=ClimMean(iidx,1)+ClimIDX(it,iidx)
         Enddo
         ClimMean(iidx,1)=ClimMean(iidx,1)/Real(Ntime)
         Do it=1, Ntime
         !   ClimIDX(it,iidx)=ClimIDX(it,iidx)-ClimMean(iidx,1)
            ClimIDX(it,iidx)=ClimIDX(it,iidx)-tmp(iidx)         ! Anomaly
         Enddo
      Enddo
      xlat(:)=1
      Lon(:)=xlon(:); Lat(:)=xlat(:)
      it1=0
      Do iyr=NYr1, NYr2
         Do it=1, NYrs
            if (Int(XYrs(it))==iyr) then
               Do imon=1, NMons
                  it1=it1+1; xtime(it1)=iyr*100+imon
               Enddo
            endif
         Enddo
      Enddo

      print*, Ntime, Nlon, Nlat
      End Subroutine Read_Clim_IDX

!------------------------------------------------------------------------------------
      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(:,:,:)
         !PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:), Mask(:,:)
         !xlon(:), xlat(:), xtime(:), Nlon, Nlat, Ntime

         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, PDSI, "SPI_03")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPEI12, "SPI_06")
         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 (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)

         print*,"In Read_Drough_IDX: Reading data ... done !"
         ! Calculate Anomaly:
         Call Calc_Anomaly_Dr (PDSI)
         Call Calc_Anomaly_Dr (SPEI12)
         Call Calc_Anomaly_Dr (SPI12)
         print*,"In Read_Drough_IDX: Calculate Anomaly ... done !"
         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) ist=ist+1
            Enddo
         Enddo
         Nst=ist
         Allocate (Dr_IDX(Nst,Ntime,NDr_idx))
         ist=0
         Do ilat=1, Nlat
            Do ilon=1, Nlon
               if (Mask(ilon,ilat)>0) then
                  ist=ist+1
                  Dr_IDX(ist,1:Ntime,1)=PDSI  (ilon,ilat,1:Ntime)
                  Dr_IDX(ist,1:Ntime,2)=SPI12 (ilon,ilat,1:Ntime)
                  Dr_IDX(ist,1:Ntime,3)=SPEI12(ilon,ilat,1:Ntime)
               endif
            Enddo
         Enddo
      End Subroutine Read_Drough_IDX
!------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------
      Subroutine Calc_Anomaly_Dr (PDSI)
         Integer :: ilon, ilat, it, Nt, imon, imon1, iyr
         Integer :: Nyear, Nmon
         Real :: tmp, xmean(12), ic(12)
         Real :: PDSI(Nlon,Nlat,Ntime)

         Nyear=30 ! 1976-2005
         Nmon =12

         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Mask(ilon,ilat)>0) then
                  Do imon=1, Nmon
                     xmean(imon)=0.; ic(imon)=0
                     Do it=1, Ntime
                        imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                        if (imon1==imon.AND.PDSI(ilon,ilat,it)/=Miss) then
                           xmean(imon)=xmean(imon)+PDSI(ilon,ilat,it); ic(imon)=ic(imon)+1
                        endif
                     Enddo
                     if (ic(imon)>0) then; xmean(imon)=xmean(imon)/ic(imon); else; xmean(imon)=Miss; endif
                  Enddo
                  Do it=1, Ntime
                     imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                     if (PDSI(ilon,ilat,it)/=Miss) PDSI(ilon,ilat,it)=PDSI(ilon,ilat,it)-xmean(imon1)
                  Enddo
               endif
            Enddo
         Enddo
      End Subroutine Calc_Anomaly_Dr
!------------------------------------------------------------------------------------
      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, Allocatable :: PCA(:,:,:), xmode(:), EOF(:,:), X(:)

         !-------------------------------------------------------------------------         
         if (.not.Allocated(PCA)) Allocate (PCA(Nlon,Nlat,nmodes))
         if (.not.Allocated(xmode)) Allocate(xmode(nmodes))
         PCA(:,:,:)=Miss; k=0
         Do j=1, Nlat
            Do i=1, Nlon
               !if (Mask(i,j)>0) then
                  k=k+1; PCA(i,j,1:nmodes) = eigenvectors(k,1:nmodes)
               !endif
            Enddo
         Enddo
         Do i=1, nmodes; 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,nmodes, 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, &
                 nmodes, trim(var_ln), "-", "lon", "lat", "mode", Miss)
         endif
         Deallocate (PCA)
         !-------------------------------------------------------------------------
         if (.not.Allocated(X)) Allocate (X(nmodes))
         X(1:nmodes) = eigenvalues(1:nmodes)
         Call wr_ncdf_1d_apd (F_pca, X, trim(IDX_nm(iidx))//"_eig", nmodes, &
                "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 )
         X(1:nmodes) = cumvariance(1:nmodes)
         Call wr_ncdf_1d_apd (F_pca, X, trim(IDX_nm(iidx))//"_cum",nmodes, &
                "Cumuluted variances of "//trim(IDX_nm(iidx)), &
                "-", "mode", Miss )
         Deallocate (X)
         !-------------------------------------------------------------------------
         if (.not.Allocated(EOF)) Allocate (EOF(nmodes,Nt))
         Do i=1, nmodes; EOF(i,1:Nt)=princomp(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, nmodes,Nt, EOF, "eof_"//trim(IDX_nm(iidx)), &
                     trim(var_ln), "-", "mode", "mon", Miss)
         else
            Call wr_ncdf_2d_apd (F_eof,EOF,"eof_"//trim(IDX_nm(iidx)), nmodes,Nt, &
                     trim(var_ln), "-", "mode", "mon", Miss)
         endif

         Deallocate (EOF)
      End Subroutine Write_EOF_PCA
!----------------------------------------------------------------------------------!
      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/Sub_Centralize.f90"
!------------------------------------------------------------------------------------
      END