Use NetCDF
      Implicit None


      Real, Parameter :: Miss=-99.0
      integer, parameter :: nmodes = 3
      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
!------------------------------------------------------------------------------------
      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="/home/tanpv/Progs/Dr_SEA/Results/Dr_PDSI_SPEI_SPI.nc"
      Call Read_Drough_IDX (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-1
         Data(1:Nx,1:Nt) = Dr_IDX(1:Nst,1:Ntime,iidx)
         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 )
      !------------------------------------------------------------------------------
          ! Write results to files:
          Fname="pca_drou_idx.nc"; Fname1="eof_timeseries_drou_idx.nc"
          Call Write_EOF_PCA (Fname,Fname1,iidx)
      Enddo


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

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