Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Yr1=1960, Yr2=2019 
      Integer, Parameter :: Nmon=12, Nreg=4+1

      ! D0 (abnormally dry, 20Ð30% percentile), D1 (moderate drought, 10Ð20%), D2 (severe drought, 5Ð10%), 
      ! D3 (extreme drought, 2Ð5%) and D4 (exceptional drought, <2%)

      Real, Parameter :: pd0=0.30, pd1=0.20, pd2=0.10, pd3=0.05, pd4=0.02   ! Probability to calculate percentiles
      Integer, parameter :: NGeo_bin=14, Nidx_bin=14
      Real, Parameter :: Geo_bin(0:NGeo_bin)=(/ 0, 2, 4, 6, 10, 14, 18, 24, 30, 36, 42, 50, 60, 76, 101/)
      Real, Parameter :: XGeo_bin(NGeo_bin) =(/   1, 3, 5, 8, 12, 16, 21, 27,  33, 39, 46, 55, 68, 80/)
      Real, Parameter :: Idx_bin(0:Nidx_bin)=(/  -1000.0,  -3.0,  -2.5,  -2.0,  -1.5,  -1.0,  -0.5,  0.0,  0.5,  1.0,  1.5,  2.0,  2.5,  3.0, 1000.0/) 
      Real, Parameter :: Xidx_bin(Nidx_bin) =     (/ -3.25,  -2.75, -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75,  3.25/)

      ! Mo ta thong tin SL dau vao:
      !---------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:), Regs(:,:), lonreg(:), latreg(:)
      Integer :: Nlon, Nlat, Ntime, Nlonreg, Nlatreg
      Logical :: Status

      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:)  ! (Nlon,Nlat,Ntime)
      Real, allocatable :: Mask(:,:) !, Mask_all(:,:)                  ! (Nlon,Nlat)
      Real, allocatable :: Dr_IDX(:,:,:)                             ! (Nreg,Ntime,NDr_idx)
      Real, allocatable :: Qp(:,:)                                   ! (Nreg,NDr_idx)
      !---------------------------------------------------------------------------------
      Real, allocatable :: Events(:,:,:)                           ! (Nreg,Ntime,NDr_idx): =1 if Dr_IDX < Qp

      Integer :: Nst, iidx, Nyear
      Integer, Parameter :: NDr_idx=2 ! SPEI12, SPI12
      Integer, Parameter :: Nsev=3    ! Moderate, Severe, Extreme
      Character (Len=3), parameter :: Sevnm(Nsev) = (/"Mod","Sev","Ext"/)
      Real :: Thres(NDr_idx,Nsev), Fix_Thres(NDr_idx) 

      !Real, parameter :: Fix_Thres(NDr_idx) = (/-2.0, -1.0, -1.0/)    ! Moderate
      !Real, parameter :: Fix_Thres(NDr_idx) = (/-3.0, -1.5, -1.5/)    ! Severe
      !Real, parameter :: Fix_Thres(NDr_idx) = (/-4.0, -2.0, -2.0/)    ! Extreme

      Logical, parameter :: Fix_Ths(NDr_idx) = (/.true.,.true./)     ! Use Fix_Thres if Fix_Ths=.true.

      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"spe12","spi"/)
                                   !(/"PDSI","SPI12","SPE12"/)
      Integer, Parameter :: NPeriod=3, NSce=2, Nmodels=8 ! 8=6-models+Ens+Obs
      Character (Len=250) :: Fout, Fname, LandMaskFile, RegionFile, Fn_mon
      Character (Len=250) :: Data_Dir, Model_nm(Nmodels), Model_Obs, Out_Dir
      !Integer :: YSta(NPeriod), YEnd(NPeriod) !, itime1(NPeriod), itime2(NPeriod)
      !Character (Len=5) :: Sce_nm(NSce), Period(NPeriod)
      !---------------------------------------------------------------------------------
      !Namelist /Infor/ Data_Dir, Period, YSta, YEnd, Model_nm, Sce_nm, LandMaskFile, RegionFile, Out_Dir
      Integer :: ilon,ilat,it,imon,iper,imod,isce
      Integer :: Nt0, itime1, itime2    ! Timeseries 
      !---------------------------------------------------------------------------------

      ! For each period: 1962-1972, 1973-1983, 1984-1994, 1995-2005, 2006-2015, 1962-2015 (6 periods)
      ! Drought frequency (DF): the number of events (spells)
      Integer, parameter :: Ndec=7
      Real, allocatable :: DF(:,:,:)       ! (nlon,nlat,ndec,ndr_idx)
      Real, allocatable :: DF_check(:,:,:)
      Real, allocatable :: DFY(:,:,:)      ! (nlon,nlat,nyear,ndr_idx)
      Real, allocatable :: DFM(:,:,:)      ! (nlon,nlat,Nmon,ndr_idx)

      ! Total drought duration (TDD): the sum of the durations
      Real, allocatable :: TDD(:,:,:)      ! (nlon,nlat,ndec,ndr_idx)
      !Real, allocatable :: TDDY(:,:,:)      ! (nlon,nlat,nyear,ndr_idx)
      ! Total drought severity (TDS): the sum of the severities of the drought events
      Real, allocatable :: TDS(:,:,:)      ! (nlon,nlat,ndec,ndr_idx)
      Real, allocatable :: TDSY(:,:,:)      ! (nlon,nlat,nyear,ndr_idx)
      ! Spatial Extent (1962-2005):
      Real, allocatable :: Geo(:,:)          ! (Ntime,ndr_idx)
      Real, allocatable :: GeoReg(:,:,:)     ! (Ntime,Nreg,ndr_idx)
      ! Areal mean of drought index:
      Real, allocatable :: Areal_mean(:,:)   ! (Ntime,ndr_idx)
      !---------------------------------------------------------
      !Trend:
      Real, allocatable :: Sen_DFY(:,:,:),  prob_DFY(:,:,:)   !(nlon,nlat,ndr_idx)
      Real, allocatable :: Sen_TDSY(:,:,:), prob_TDSY(:,:,:)
      Real, allocatable :: Sen_GeoReg(:,:), prob_GeoReg(:,:)  ! (Nreg,ndr_idx)
      !---------------------------------------------------------
      Real, allocatable :: Freq_idx(:,:,:)     !(Nidx_bin,Ndec,NDr_idx) idx frequency
      Integer :: Nsec, Nyrsec, isev, Nsec1

      Thres(1:NDr_idx,1) = (/-1.0, -1.0/)
      Thres(1:NDr_idx,2) = (/-1.5, -1.5/) 
      Thres(1:NDr_idx,3) = (/-2.0, -2.0/)

      !LandMaskFile="/work/users/hmo/tanpv/CORDEX_SEACLID/Land_mask_sea.nc"
      !Call Read_Land_Mask

      !RegionFile="/home/tanpv/Progs/Dr_SEA/Post_APH/NC1/kmeans_4.nc"
      !Call rd_ncdf_2d (RegionFile,lonreg,latreg,Nlonreg,Nlatreg,Regs, "cluster")

      Nyear = Yr2-Yr1+1           ! Bo nam dau tien neu Yr1=1959
      Nyrsec = Nyear/(Ndec-1)   ! Moi giai doan (20 nam)
      
      !Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_APH_SEA_61-15.nc" !BL_APH_SEA.nc"
      Fname="/work/projects/G09/SEACLID/GPCC_CRU_SEA/SPEI_SPI_GPCC_CRU_Reg_59_19_New.nc"  !/BL_APH_SEA_Reg_61-15.nc"
      Call Read_Drough_IDX (Fname)    ! Return Dr_IDX(nlon,nlat,Ntime,NDr_idx), xlon, xlat

      Allocate (DF_check (Nreg,ndec,ndr_idx)) ! Check if DF is suitable

      Do isev=1,Nsev
         Fix_Thres(1:NDr_idx) = Thres(1:NDr_idx,isev)

         Call Calc_Per_Events (pd3)      ! if Fix_Ths=.false. then used pd3=5% --> Events(ilon,ilat,it,iidx)
         Nsec = (Ntime)/(Ndec-1)    ! 11 year = 132 mon
         Nsec1 = Nsec !-Nmon               ! 10 year

         Call Calc_chars
         !Call Calc_Trend_DFY_TDSY_GeoReg
         if (isev>1) Where (DF/=Miss.AND.DF_check/=Miss.AND.DF>DF_check) DF=DF_check

         ! Write out:
         Fname="/home/tanpv/Progs/2021/NC/CRU_Chars_Reg_"//Sevnm(isev)//"_10.nc"      ! Moderate
         Fn_mon="/home/tanpv/Progs/2021/NC/CRU_Fre_mon_Reg_"//Sevnm(isev)//"_10.nc"
         Fout="/home/tanpv/Progs/2021/NC/CRU_Chars_Yearly_Reg_"//Sevnm(isev)//"_10.nc"

         Call Write_Char_to_NC (Fname, Fn_mon, Fout)

         DF_check = DF   
         Deallocate (Qp,Events)
         Deallocate (DF, TDD, TDS, DFY, TDSY, DFM) !, Geo, GeoReg, Freq_idx, Areal_mean)
         !Deallocate (Sen_DFY,prob_DFY,Sen_TDSY,prob_TDSY,Sen_GeoReg,prob_GeoReg)
      Enddo
   CONTAINS
   !---------------------------------------------------------------------------------
   Subroutine Calc_chars
      Integer :: ilon,ilat,it,iidx
      Real :: X(Ntime), X1(Nsec), X2(Nmon), M1(Ntime), IMAX1, Spells(Ntime), N_Spell_1, I_Start(Ntime), I_End(Ntime)
      Integer :: ireg, N, it1,it2, imon, imon1, i, k, N_Spell, Imax, ista, iend, Nsec0, Nyrsec0
      Real :: tmp, ntmp,tmp1, ntmp1, Smean, Smax, In_mean, In_max, Amean, Se0, Sem, Sex, Inm, Inx

      Allocate (DF (Nreg,ndec,ndr_idx))
      Allocate (TDD(Nreg,ndec,ndr_idx))
      Allocate (TDS(Nreg,ndec,ndr_idx))

      Allocate (DFY (Nreg,Nyear,ndr_idx))
      !Allocate (TDDY(nlon,nlat,Nyear,ndr_idx))
      Allocate (TDSY(Nreg,Nyear,ndr_idx))
      Allocate (DFM(Nreg,Nmon,ndr_idx))

      !Allocate (Geo(Ntime,ndr_idx))
      !Allocate (GeoReg(Ntime,Nreg,ndr_idx))
      !Allocate (Freq_idx(Nidx_bin,ndec,NDr_idx))
      !Allocate (Areal_mean(Ntime,ndr_idx))

      DF=Miss; TDD=Miss; TDS=Miss; !Geo=Miss; GeoReg=Miss
      M1=Miss; Spells=Miss; I_Start=Miss;I_End=Miss
      !Areal_mean=Miss; 
      DFY=Miss; TDSY=Miss; DFM=Miss

      Do iidx=1,NDr_idx
         Do ireg=1,Nreg !ilon=1, Nlon
            !Do ilat=1, Nlat
               !if (Mask(ilon,ilat) /= Miss) then
                  X(:) = Events(ireg,:,iidx)
                  Call NSpell_and_Leng (X, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss)
                  N_Spell=NInt(N_Spell_1)
                  DF (ireg,ndec,iidx) = N_Spell_1/Nyear*10.  ! per decade               ! Annual average
                  TDD(ireg,ndec,iidx) = Sum(Spells(1:N_Spell))/Nyear*10.  ! per decade 
                  tmp=0
                  Do it=1,Ntime
                     if (X(it)==1.AND.Dr_IDX(ireg,it,iidx)/=Miss) tmp = tmp + Dr_IDX(ireg,it,iidx)
                  Enddo
                  TDS(ireg,ndec,iidx) = tmp/Nyear*10.  ! per decade

                  Do k=1,ndec-1
                     it1=(k-1)*nsec+1
                     !if (k<ndec-1) then
                        it2=it1+nsec-1; Nsec0=Nsec; Nyrsec0=Nyrsec
                     print*,it1,it2,nsec,Nyrsec
                     !else
                     !   it2=it1+nsec1-1; Nsec0=Nsec1; Nyrsec0=Nyrsec-1
                     !endif
                     X1(1:Nsec0) = Events(ireg,it1:it2,iidx)
                     Call NSpell_and_Leng (X1, Nsec0, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss)
                     N_Spell=NInt(N_Spell_1)
                     DF (ireg,k,iidx) = N_Spell_1/Nyrsec0*10.  ! per decade               ! Annual average
                     TDD(ireg,k,iidx) = Sum(Spells(1:N_Spell))/Nyrsec0*10.  ! per decade
                     tmp=0; i=0
                     Do it=it1,it2   ! So thang trong thap ky
                        i=i+1
                        if (X1(i)==1.AND.Dr_IDX(ireg,it,iidx)/=Miss) tmp = tmp + Dr_IDX(ireg,it,iidx)
                     Enddo
                     TDS(ireg,k,iidx) = tmp/Nyrsec0*10.  ! per decade
                  Enddo
                  ! Yearly:
                  Do iyr=1,Nyear
                     it1=(iyr-1)*Nmon+1; it2=it1+Nmon-1
                     !print*,it1,it2
                     X2(1:Nmon) = Events(ireg,it1:it2,iidx)
                     tmp=0; tmp1=0; i=0
                     Do it=it1,it2   ! So thang trong nam iyr
                        i=i+1
                        if (X2(i)==1.AND.Dr_IDX(ireg,it,iidx)/=Miss) then
                           tmp=tmp+1; tmp1=tmp1+Dr_IDX(ireg,it,iidx)      ! So thang han va tong Dr_IDX<Thres
                        endif
                     Enddo
                     DFY (ireg,iyr,iidx) = tmp; TDSY(ireg,iyr,iidx) = tmp1
                     !print*,DFY (ilon,ilat,iyr,iidx),TDSY(ilon,ilat,iyr,iidx)
                  Enddo 
                  !print"(10F7.3)",DFY (ilon,ilat,1:Nyear,iidx)
               !endif
            !Enddo
         Enddo
      Enddo 
      ! Monthly Freq (%):
      Do iidx=1,NDr_idx
         Do ireg=1,Nreg !ilon=1, Nlon
            !Do ilat=1, Nlat
            !   if (Mask(ilon,ilat) /= Miss) then
                  Do imon=1,Nmon
                     ntmp=0
                     Do it=1,Ntime
                        if (Events(ireg,it,iidx)==1) then
                           imon1=Mod(it,Nmon); if (imon1==0) imon1=Nmon
                           if (imon1==imon) ntmp=ntmp+1   
                        endif
                     Enddo
                     DFM(ireg,imon,iidx)=real(ntmp)/real(Nyear)*100
                  Enddo
            !   endif
            !Enddo
         Enddo
      Enddo

   End Subroutine Calc_chars
   !---------------------------------------------------------------------------------
   Subroutine Write_Char_to_NC(Fname, Fn_fr_mon, Fn_yrs)
      Character (Len=*) :: Fname, Fn_fr_mon, Fn_yrs
      Character (Len=50)  :: var_n, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=250) :: var_ln
      Real :: mons(Nmon), x_idx(NDr_idx), xmon(Ntime), xreg(Nreg), xdec(Ndec), xreg1(Nreg+1)
      Real :: year(Nyear)
      Real, Allocatable, Dimension (:,:,:) :: Xtmp
      Integer :: imon, ireg, iidx, iyr,idec

!      Allocate (TDS(Nreg,ndec,ndr_idx))

!      Allocate (DFY (Nreg,Nyear,ndr_idx))
      !Allocate (TDDY(nlon,nlat,Nyear,ndr_idx))
!      Allocate (TDSY(Nreg,Nyear,ndr_idx))
!      Allocate (DFM(Nreg,Nmon,ndr_idx))

      lon_nm="reg"; lat_nm="decade"; t_nm="idx";
      t_unit="years since 1960-01-01 "

      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Do iidx=1,NDr_idx; x_idx(iidx)=iidx; Enddo
      Do imon=1,Ndec; xdec(imon)=imon; Enddo 
      Do imon=1,Nyear; year(imon)=Yr1+imon-1; Enddo
      !xdec(1)=1960; xdec(2)=1980; xdec(3)=2000; xdec(4)=2019; ! xdec(5)=2005; xdec(6)=2015

      ! Drought duration:
      var_ln="Annual frequency (DF) for 62-72,73-83,84-94,95-05,06-15 and 1962-2015. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="spell"
      Call wr_ncdf_3d (Fname,xreg,xdec,x_idx,Nreg,Ndec,NDr_idx,DF,"df",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      print*,"Ndec=",Ndec
      Do imon=1,Ndec; print*,"idec=",imon, DF(:,imon,1); Enddo

      var_ln="Annual drought duration (TDD) for 62-72,73-83,84-94,95-05,06-15 and 1962-2015. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="mon"
      Call wr_ncdf_3d_apd (Fname,TDD,"tdd",Nreg,Ndec,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      var_ln="Annual drought drought severity (TDS) for 62-72,73-83,84-94,95-05,06-15 and 1962-2015. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="-"
      Call wr_ncdf_3d_apd (Fname,TDS,"tds",Nreg,Ndec,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Monthly Frequency:
      var_ln="Monthly frequency of drought (dfm) for  period 1962-2015. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="%"; lat_nm="month"
      Do imon=1,Nmon; mons(imon)=imon; Enddo
      Call wr_ncdf_3d (Fn_fr_mon,xreg,mons,x_idx,Nreg,Nmon,NDr_idx,DFM,"dfm",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Yearly:
      var_ln="Yearly number of drought months"
      var_unit="mon"; lat_nm="year";
      Call wr_ncdf_3d (Fn_yrs,xreg,year,x_idx,Nreg,Nyear,NDr_idx,DFY,"dfy",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Yearly drought severity"; var_unit="1"
      Call wr_ncdf_3d_apd (Fn_yrs,TDSY,"tdsy",Nreg,Nyear,NDr_idx,var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Write to TXT:
      Open (3,file="/home/tanpv/Progs/2021/NC/Char_Reg_"//Sevnm(isev)//".txt",status="Unknown")
      Do iidx=1,NDr_idx
         write (3,*) "Drought IDX: ",IDX_nm(iidx)
         write (3,"(8A8)") "Period","60-69","70-79","80-89","90-99","00-09","10-19","60-19"
         Do ireg=1,Nreg
            write (3,"(A,I1)") "Sub-region R",ireg
            write (3,"(A8,7F8.3)") "DF",  ( DF(ireg,idec,iidx),idec=1,Ndec)
            write (3,"(A8,7F8.3)") "TDD", (TDD(ireg,idec,iidx),idec=1,Ndec)
            write (3,"(A8,7F8.3)") "TDS", (TDS(ireg,idec,iidx),idec=1,Ndec)
         Enddo
         write (3,*)
         write (3,"(A8,12(2x,A4,I2.2))") "Mon",("Mon_",imon,imon=1,Nmon)
         Do ireg=1,Nreg
            write (3,"(6x,A,I1,12F8.2)") "R",ireg, (DFM(ireg,imon,iidx),imon=1,Nmon)
         Enddo
      Enddo
      Close (1)
      
   End Subroutine Write_Char_to_NC
   !---------------------------------------------------------------------------------
   Subroutine Calc_Per_Events (p)
      ! Return: Qp and Events
      Integer :: ilon, ilat, iidx, it, ireg
      Real :: p
      Real :: X(Ntime)
         
      Allocate (Qp(Nreg,NDr_idx),Events(Nreg,Ntime,NDr_idx))

      Events = Miss; Qp = Miss   ! Need to be check !!!

      Do iidx=1,NDr_idx
         Do ireg=1,Nreg !ilon=1, Nlon
            !Do ilat=1, Nlat 
            !   if (Mask(ilon,ilat) /= Miss) then
                  if (Fix_Ths(iidx)) then                ! Use Fix_Thres if Fix_Thrs=.True.
                     Qp(ireg,iidx) = Fix_Thres(iidx)
                  else
                     X(:) = Dr_IDX(ireg,:,iidx)
                     Qp(ireg,iidx) = percentile (X,Ntime,p,Miss)
                  endif
                  Do it=1, Ntime
                     if (Dr_IDX(ireg,it,iidx) /=Miss .AND. Dr_IDX(ireg,it,iidx) <= Qp(ireg,iidx)) then 
                        Events(ireg,it,iidx) = 1
                     !else if (Dr_IDX(ireg,it,iidx)==Miss) then
                     !   Events(ireg,it,iidx) = Miss
                     else if (Dr_IDX(ireg,it,iidx) /=Miss .AND. Dr_IDX(ireg,it,iidx)>Qp(ireg,iidx)) then
                        Events(ireg,it,iidx) = 0
                     endif               
                  Enddo
            !   endif
            !Enddo
         Enddo
      Enddo

   End Subroutine Calc_Per_Events
   !---------------------------------------------------------------------------------
      Subroutine Read_Drough_IDX (Fname)
         Integer :: ist, ilon, ilat, it, it1, it2, Nt, imon, imon1, iyr
         Integer :: nlat1, nlat0, nlon, nlat
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         Real, Allocatable :: X0(:), X1(:,:), X2(:,:), X3(:,:)
         Real, Allocatable :: xlon(:), xlat(:)

         !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, "spe12")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "spi")
         Deallocate (xlon, xlat, xtime)

         !Where (IsNaN(PDSI)); PDSI=Miss; Endwhere
         Where (IsNaN(SPEI12)); SPEI12=Miss; Endwhere
         Where (IsNaN(SPI12)); SPI12=Miss; Endwhere

         Ntime = Nt-Nmon     ! Skip first year due to SPI and SPEI were missed
         it1 = Nmon+1; it2 = Nt
         print*,"DrIDX: ",Nt,Ntime

         !Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime))
         Allocate (X1(Nreg,Ntime),X2(Nreg,Ntime),X3(Nreg,Ntime))
         !X1(:,:,1:Ntime)=PDSI  (:,:,it1:it2) !Nmon+1:Nt)
         X2(:,1:Ntime)=SPEI12(:,1,it1:it2) !Nmon+1:Nt)
         X3(:,1:Ntime)=SPI12 (:,1,it1:it2) !Nmon+1:Nt)

         Deallocate (SPEI12,SPI12)
         !Allocate (PDSI(nlon,nlat,Ntime),SPEI12(nlon,nlat,Ntime),SPI12(nlon,nlat,Ntime))

         !PDSI=X1; SPEI12=X2; SPI12=X3

         Allocate (Dr_IDX(Nreg,Ntime,NDr_idx))

         Dr_IDX(:,:,:) = Miss
         Do ist=1,Nreg
            if (.not.Empty_Ar1(X2(ist,1:Ntime),Ntime).AND. &
                .not.Empty_Ar1(X3(ist,1:Ntime),Ntime)) then
               Dr_IDX(ist,:,1)=X2(ist,:)
               Dr_IDX(ist,:,2)=X3(ist,:)
            endif
         Enddo

         Deallocate (X1, X2, X3)
         ! Note: Nreg=Nlon,

         !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(ilon,ilat,:,1)=PDSI(ilon,ilat,:)
         !            Dr_IDX(ilon,ilat,:,2)=SPI12(ilon,ilat,:)
         !            Dr_IDX(ilon,ilat,:,3)=SPEI12(ilon,ilat,:)
         !         else
         !            Mask(ilon,ilat)=Miss
         !         endif
         !      endif
         !   Enddo
         !Enddo
         !Nst=ist   ! Number of gridpoints
         !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 "/home/tanpv/Mylibs/Sub_Percentile.f90"
   Include "/home/tanpv/Mylibs/Sub_NSpell_and_Leng.f90"
   Include "/home/tanpv/Mylibs/Ave.f90"
   Include "/home/tanpv/Mylibs/Stdev.f90"
   Include "/home/tanpv/Mylibs/Sub_Trend_MK_Linear.f90"
   !Include "/home/tanpv/Mylibs/Sum.f90"

   END Program