Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Yr1=1976, Yr2=2005, Nmon=12

      ! 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

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

      Real, allocatable :: PDSI(:,:,:), SPI12(:,:,:), SPEI12(:,:,:)  ! (Nlon,Nlat,Ntime)
      Real, allocatable :: Mask(:,:), Mask_all(:,:)                  ! (Nlon,Nlat)
      Real, allocatable :: Dr_IDX(:,:,:,:)                           ! (Nlon,Nlat,Ntime,NDr_idx)
      Real, allocatable :: Qp(:,:,:)                                 ! (Nlon,Nlat,NDr_idx)
      !---------------------------------------------------------------------------------
      Real, allocatable :: Events(:,:,:,:)                           ! (Nlon,Nlat,Ntime,NDr_idx): =1 if Dr_IDX < Qp
      Real, allocatable :: FrMo(:,:,:,:)                             ! (Nlon,Nlat,Nmon,NDr_idx)-% of Event/mon
      Real, allocatable :: Freq(:,:,:,:)                             ! (Nlon,Nlat,Nyear,NDr_idx)-Num of Event/Year
      Real, allocatable :: Geo(:,:)                                  ! (Ntime,NDr_idx) - Geo Ext in each mon
      Real, allocatable :: Dur_M(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Mean Duration
      Real, allocatable :: Dur_X(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Max Duration
      Real, allocatable :: Sev_M(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Mean Severity
      Real, allocatable :: Sev_X(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Max Severity
      Real, allocatable :: Inten_M(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Mean Intensity
      Real, allocatable :: Inten_X(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Max Intensity
      Real, allocatable :: Inter_M(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Mean Inter-arrival
      Real, allocatable :: Inter_X(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Max Inter-arrival

      Real, allocatable :: Fr_Max(:,:,:)                           ! (Nlon,Nlat,NDr_idx) - Maximum annual Freq (mon/year) 
      Real, allocatable :: Fr_Std(:,:,:)                           ! (Nlon,Nlat,NDr_idx) - Sdt annual Freq (mon/year)
      !Real, allocatable :: q25(:,:,:)                             ! (Nlon,Nlat,NDr_idx) - Freq q25 (mon/year)
      !Real, allocatable :: q75(:,:,:)                             ! (Nlon,Nlat,NDr_idx) - Freq q75 (mon/year)
      !Real, allocatable :: IQR(:,:,:)                             ! (Nlon,Nlat,NDr_idx) - Freq IQR (mon/year)

      Integer :: Nst, iidx, Nyear
      Integer, Parameter :: NDr_idx=3 ! PDSI, SPEI12, SPI12
      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"pdsi","spi12","spe12"/)
                                   !(/"PDSI","SPI12","SPE12"/)
      Character (Len=100) :: Fname
      !---------------------------------------------------------------------------------

   Fname="Data/BL_CNRM-CM5.nc" ! "Dr_PDSI_SPEI_SPI.nc"
   Call Read_Drough_IDX (Fname) 
   Nyear = Ntime/Nmon
   Allocate (Qp(Nlon,Nlat,NDr_idx), Events(Nlon,Nlat,Ntime,NDr_idx))
   Events (:,:,:,:) = Miss
   Call Calc_Per_Events (pd3)  ! 0.05

   Allocate (Freq(Nlon,Nlat,Nyear,NDr_idx),FrMo(Nlon,Nlat,Nmon,NDr_idx))
   Allocate (Geo(Ntime,NDr_idx))
   Allocate (Dur_M(Nlon,Nlat,NDr_idx), Dur_X(Nlon,Nlat,NDr_idx))
   Allocate (Sev_M(Nlon,Nlat,NDr_idx), Sev_X(Nlon,Nlat,NDr_idx))
   Allocate (Inten_M(Nlon,Nlat,NDr_idx), Inten_X(Nlon,Nlat,NDr_idx))
   Allocate (Inter_M(Nlon,Nlat,NDr_idx), Inter_X(Nlon,Nlat,NDr_idx))

   !Allocate (Mean(Nlon,Nlat,NDr_idx), q25(Nlon,Nlat,NDr_idx), q75(Nlon,Nlat,NDr_idx), IQR(Nlon,Nlat,NDr_idx))
   Allocate (Fr_max(Nlon,Nlat,NDr_idx), Fr_Std(Nlon,Nlat,NDr_idx))

   ! Set all to Miss:
   Freq=Miss; FrMo=Miss; Geo=Miss; Dur_M=Miss; Dur_X=Miss
   Sev_M=Miss; Sev_X=Miss; Inten_M=Miss; Inten_X=Miss; Inter_M=Miss; Inter_X=Miss
   Fr_Max=Miss; Fr_Std=Miss !q25=Miss; q75=Miss; IQR=Miss

   Call Calc_Chars
   Call Write_Fre_to_NC ("Freq.nc","Fre_mon.nc")
   Call Write_Char_to_NC("Chars.nc","Geo.nc")

   CONTAINS
   !---------------------------------------------------------------------------------
   Subroutine Write_Char_to_NC(Fname,Fn_geo)
      Character (Len=*) :: Fname, Fn_geo
      Character (Len=50) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
      Real :: x_idx(NDr_idx), xmon(Ntime)

      lon_nm="lon"; lat_nm="lat"; t_nm="idx";
      t_unit="years since 1977-01-01 "
      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3
      ! Drought duration:
      var_ln="Mean duration of drought (month). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="mon"
      Call wr_ncdf_3d (Fname,xlon,xlat,x_idx,nlon,nlat,NDr_idx, Dur_M, "dur_m",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Maximum duration of drought (month). IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d_apd (Fname, Dur_X,"dur_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      ! Drought severity:
      var_ln="Mean severity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="-"
      Call wr_ncdf_3d_apd (Fname, Sev_M,"sev_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Maximum severity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d_apd (Fname, Sev_X,"sev_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      ! Drought intensity:
      var_ln="Mean intensity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d_apd (Fname, Inten_M,"inten_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Maximum intensity of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d_apd (Fname, Inten_X,"inten_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      ! Drought inter-arrival:
      var_ln="Mean inter-arrival of drought"; var_unit="mon"
      Call wr_ncdf_3d_apd (Fname, Inter_M,"inter_m",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Maximum inter-arrival of drought. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d_apd (Fname, Inter_X,"inter_x",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Drought characteristics' frequency (Mean, q25, q75, IQR of Freq):
      var_ln="Maximum Annual drought frequency (mon/year)"
      var_unit="mon/year"
      Call wr_ncdf_3d_apd (Fname, Fr_Max,"fr_max",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Standard deviation of Annual drought frequency (mon/year)"
      Call wr_ncdf_3d_apd (Fname, Fr_Std,"fr_std",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      !var_ln="First quartile (q25) of annual drought frequency (mon/year)"
      !Call wr_ncdf_3d_apd (Fname, q25,"fr_q25",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      !var_ln="Third quartile (q75) of annual drought frequency (mon/year)"
      !Call wr_ncdf_3d_apd (Fname, q75,"fr_q75",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      !var_ln="Interquartile (IQR) of annual drought frequency (mon/year)"
      !Call wr_ncdf_3d_apd (Fname, IQR,"fr_iqr",nlon,nlat,NDr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Geographic Extent:
      Do iyr=1, Ntime; xmon(iyr)=iyr; Enddo
      var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1977. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="%"; lon_nm="mon"; lat_nm="idx"
      Call wr_ncdf_2d (Fn_geo,xmon,x_idx,Ntime,NDr_idx, Geo,"geo", var_ln, var_unit, lon_nm, lat_nm, Miss)

   End Subroutine Write_Char_to_NC
   !---------------------------------------------------------------------------------
   Subroutine Write_Fre_to_NC (Fname1, Fname2)
      Character (Len=*) :: Fname1, Fname2
      Character (Len=50) :: var_n, var_ln, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm 
      Real :: x_idx(NDr_idx), xyear(Nyear), xmon(Nmon)
      Integer :: iyr, iidx
      Real, Allocatable :: X(:,:,:,:) !Nlon,Nlat,NDr_idx,Nyear)
      ! Freq:
      lon_nm="lon"; lat_nm="lat"; t_nm="year"; z_nm="idx"
      t_unit="years since 1977-01-01 "
      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3
      Do iyr=1, Nyear; xyear(iyr)=Yr1+iyr; Enddo
      Allocate (X(Nlon,Nlat,NDr_idx,Nyear))
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            Do iidx=1, NDr_idx
                Do iyr=1, Nyear
                   X(ilon,ilat,iidx,iyr) = Freq(ilon,ilat,iyr,iidx)
                Enddo
            Enddo
         Enddo
      Enddo
      var_ln="Number of drought month in the year. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="mon"
      print*, "Writing to file ..."
      Call wr_ncdf_4d (Fname1,xlon,xlat,x_idx,xyear,nlon,nlat,NDr_idx,Nyear, X, "freq", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
      !FrMo:
      t_nm="mon"; t_unit="months since 1977-01-01 "
      Do iyr=1, Nmon; xmon(iyr)=iyr; Enddo      ! months in the year
      Deallocate(X); Allocate (X(Nlon,Nlat,NDr_idx,Nmon))
      Do ilon=1, Nlon
         Do ilat=1, Nlat
            Do iidx=1, NDr_idx
               X(ilon,ilat,iidx,:) = FrMo(ilon,ilat,:,iidx)
            Enddo
         Enddo
      Enddo
      var_ln="Monthly Drought frequency (%): Drought frequency in each month. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="%"
      Call wr_ncdf_4d (Fname2,xlon,xlat,x_idx,xmon,nlon,nlat,NDr_idx,Nmon, X, "freq", &
                      t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
   End Subroutine Write_Fre_to_NC
   !---------------------------------------------------------------------------------
   Subroutine Calc_Chars
      Integer :: ilon, ilat, iidx, it, iyr, imon
      Real :: X(Ntime), M1(Ntime), IMAX1, Spells(Ntime), N_Spell_1, I_Start(Ntime), I_End(Ntime)
      Integer :: N, i, k, N_Spell, Imax, ista, iend
      Real :: tmp, ntmp, Smean, Smax, In_mean, In_max, Amean
      Real :: XX(Nyear)

      M1=Miss; Spells=Miss; I_Start=Miss;I_End=Miss
      Do iidx=1,NDr_idx
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Mask_all(ilon,ilat) /= Miss) then
                  X(:) = Events(ilon,ilat,:,iidx)
                  Call NSpell_and_Leng (X, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss)
                  print"(A,20I6)","iidx,Imax1,Maxval(Spells), N_Spell: ", iidx,int(Imax1),int(Maxval(Spells)), Int(N_Spell_1)
                  N_Spell=Int(N_Spell_1); Imax=Int(Imax1)  ! So dot va Do dai dot dai nhat
                  tmp=0
                  Do i=1,N_Spell
                     tmp=tmp+Spells(i)                     ! Tong do dai tat ca cac dot
                  Enddo 
                  if (N_Spell>0) then 
                     Dur_M(ilon,ilat,iidx)=tmp/N_Spell_1; Dur_X(ilon,ilat,iidx)=Imax1
                  else 
                     Dur_M(ilon,ilat,iidx)=Miss; Dur_X(ilon,ilat,iidx)=Miss
                  endif  ! Do dai trung binh & Do dai dot dai nhat
                  ! Severity and Intensity:
                  Smean=0; Smax=999.; In_mean=0; In_max=999.
                  Do i=1, N_Spell
                     ista=Int(I_Start(i)); iend=Int(I_End(i))
                     Amean=0;
                     Do k=ista, iend     ! Thoi diem bat dau/ket thuc dot thu i
                        Amean = Amean + Dr_IDX(ilon,ilat,k,iidx) ! Seve dot thu i = Tong Dr_IDX trong dot thu i
                     Enddo
                     Smean = Smean+Amean ! Tong Dr_IDX trong tat ca cac dot (Tong Seve)
                     if (Smax>Amean) Smax = Amean   ! Max Severity cua cac dot
                     Amean = Amean/Spells(i)        ! Intensity dot thu i
                     In_mean = In_mean+Amean        ! Tong Intensity cua cac dot
                     if (In_max>Amean) In_max = Amean  ! Max Intensity cua cac dot
                  Enddo
                  if (N_Spell>0) then 
                     Sev_M(ilon,ilat,iidx) = Smean/N_Spell_1; Sev_X(ilon,ilat,iidx) = Smax
                     Inten_M(ilon,ilat,iidx) = In_mean/N_Spell_1; Inten_X(ilon,ilat,iidx) = In_max
                  else
                     Sev_M(ilon,ilat,iidx) = Miss; Sev_X(ilon,ilat,iidx) = Miss
                     Inten_M(ilon,ilat,iidx) = Miss; Inten_X(ilon,ilat,iidx) = Miss
                  endif ! Seve & Inten
                  ! Inter-arrival:
                  if (N_Spell>1) then
                     Amean=0; tmp=0
                     Do i=1, N_Spell-1
                        ntmp = I_Start(i+1)-I_Start(i)   ! Khoang giua hai lan Onset lien tiep
                        Amean = Amean+ntmp               ! Tong cac Inter-Arrival
                        if (ntmp>tmp) tmp = ntmp         ! Max Inter-Arrival
                     Enddo
                     Inter_M(ilon,ilat,iidx) = Amean/(N_Spell_1-1); Inter_X(ilon,ilat,iidx) = tmp
                  else
                     Inter_M(ilon,ilat,iidx) = Miss; Inter_X(ilon,ilat,iidx) = Miss
                  endif
                  ! Freq: So events/nam (So thang han/nam)
                  it=0
                  Do iyr=1, Nyear
                     Freq(ilon,ilat,iyr,iidx) = 0
                     Do imon=1, Nmon
                        it=it+1; Freq(ilon,ilat,iyr,iidx) = Freq(ilon,ilat,iyr,iidx) + Events(ilon,ilat,it,iidx)
                     Enddo
                  Enddo

                  ! Chars of Freq: Mean, q25, q75, IQR

                  !XX(1:Nyear) = Freq(ilon,ilat,1:Nyear,iidx)
                  !X(1:Ntime) = Events(ilon,ilat,1:Ntime,iidx)
                  !print*,xlon(ilon),xlat(ilat),Int(Sum(X,mask=X>=0))
                  !Mean(ilon,ilat,iidx) = Sumx(X, Ntime, Miss)!/Nyear !Sum(X,mask=X>=0)/Nyear    ! Trung binh Mon/Year
                  !print*,Mean(ilon,ilat,iidx)

                  XX(1:Nyear) = Freq(ilon,ilat,1:Nyear,iidx)
                  Fr_Max(ilon,ilat,iidx) = Maxval(XX)         ! Nam co so thang han nhieu nhat
                  Fr_Std(ilon,ilat,iidx) = Stdev (XX,Nyear,Miss) 

                  !q25(ilon,ilat,iidx)    = Stdev (X,Ntime,Miss)      ! 
                  !q75(ilon,ilat,iidx)    = Maxval (XX(1:Nyear))
                  !if (q25(ilon,ilat,iidx)/=0) then 
                  !   IQR(ilon,ilat,iidx)=q25(ilon,ilat,iidx)/Mean(ilon,ilat,iidx)
                  !else
                  !   IQR(ilon,ilat,iidx)=Miss
                  !endif 

                  ! FrMo: Tan suat xuat hien trong tung thang (%/mon)
                  Do imon=1, Nmon
                     FrMo(ilon,ilat,imon,iidx) = 0
                     Do it=1, Ntime
                        imon1=mod(it,Nmon); if (imon1==0) imon1 = Nmon
                        if (imon1==imon) then
                           FrMo(ilon,ilat,imon,iidx) = FrMo(ilon,ilat,imon,iidx) + Events(ilon,ilat,it,iidx)
                        endif
                     Enddo
                     FrMo(ilon,ilat,imon,iidx) = FrMo(ilon,ilat,imon,iidx)/Real(Nyear)*100.0
                  Enddo
               endif
            Enddo
         Enddo
      Enddo
      ! Geo Ext:
      Do iidx=1,NDr_idx
         Do it=1, Ntime
            Geo(it,iidx) = 0
            Do ilon=1, Nlon
               Do ilat=1, Nlat
                  if (Mask_all(ilon,ilat) /= Miss) then
                     Geo(it,iidx) = Geo(it,iidx) + Events(ilon,ilat,it,iidx) 
                  endif
               Enddo
            Enddo
            Geo(it,iidx)=Geo(it,iidx)/Real(Nst)*100.0
         Enddo
         !print*,Geo(:,iidx);pause
      Enddo
   End Subroutine Calc_Chars
   !---------------------------------------------------------------------------------
   Subroutine Calc_Per_Events (p)
      ! Return: Qp and Events
      Integer :: ilon, ilat, iidx, it
      Real :: p
      Real :: X(Ntime)
         
      Do iidx=1,NDr_idx
         Do ilon=1, Nlon
            Do ilat=1, Nlat 
               if (Mask_all(ilon,ilat) /= Miss) then
                  X(:) = Dr_IDX(ilon,ilat,:,iidx)
                  Qp(ilon,ilat,iidx) = percentile (X,Ntime,p,Miss)
                  Do it=1, Ntime
                     if (Dr_IDX(ilon,ilat,it,iidx) < Qp(ilon,ilat,iidx)) then 
                        Events(ilon,ilat,it,iidx) = 1
                     else
                        Events(ilon,ilat,it,iidx) = 0
                     endif               
                  Enddo
                  !print "(30F6.2)",X(:);print*
                  !print "(30F6.2)", Qp(ilon,ilat,iidx);print*
                  !print "(30I6)", Int(Events(ilon,ilat,:,iidx)); print*
                  !X(:)=Events(ilon,ilat,:,iidx)
                  !Do it=1,Ntime; if (X(it)>0) print*,it," ",Dr_IDX(ilon,ilat,it,iidx),Qp(ilon,ilat,iidx),X(it); Enddo
                  !print*,Int(Sum(X,mask=X>0))
                  !print*,Int(sum(Events(ilon,ilat,1:Ntime,iidx))); pause

               endif
            Enddo
         Enddo
      Enddo

   End Subroutine Calc_Per_Events
   !---------------------------------------------------------------------------------
      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)
         Fname="/work/projects/G09/SEACLID/Land_mask_sea.nc"
         Call rd_ncdf_2d (trim(Fname), xlon,xlat, nlon,nlat,Mask,"land")
         Nyear=30 ! 1976-2005
         Nmon =12
         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:Ntime)=PDSI  (:,:,Nmon+1:Nt)
         X2(:,:,1:Ntime)=SPEI12(:,:,Nmon+1:Nt)
         X3(:,:,1:Ntime)=SPI12 (:,:,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
         Allocate (Dr_IDX(nlon,nlat,Ntime,NDr_idx))
         Dr_IDX(:,:,:,:) = Miss
         !Dr_IDX(:,:,:,1)=PDSI(:,:,:); Dr_IDX(:,:,:,2)=SPI12(:,:,:); Dr_IDX(:,:,:,3)=SPEI12(:,:,:)

         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
                     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_all(ilon,ilat)=Miss
                  endif
               endif
            Enddo
         Enddo
         Nst=ist   ! Number of gridpoints
         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/Sum.f90"

   END Program