Use NetCDF
      Real, Parameter :: Miss=-99.
      Integer, Parameter :: Yr1=1980, Yr2=2020 
      Integer, Parameter :: Nmon=12, Nreg=7, Nidx=24

      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/)

      Character (Len=15), Parameter :: Clim_nm(Nidx)=(/ &
                "nino12", "nino3", "nino34", "nino4", "oni", "tni", &   ! 6
                "whwp", "ammsst", "pdo", "NTA_ersst", "dmi",  &   ! 5
                "tna", "amonus", &
                "pna", "censo", "soi", "reqsoi", "np", "repac_slpa", & ! 6 
                "rindo_slpa", "qbo30", "wp", "ao", "nao"/)

      ! Mo ta thong tin SL dau vao:
      !---------------------------------------------------------------------------------
      Real, allocatable :: xlon(:), xlat(:), xtime(:), Reg(:)
      Character (Len=20), Allocatable :: Stname(:) 
      Integer :: Nst, Ntime
      Logical :: Status

      Real, allocatable :: SPI(:,:,:), SPEI(:,:,:)  ! (Nst,Ntime,Nlen)
      Real, allocatable :: Dr_IDX(:,:,:,:)                           ! (Nst,Ntime,Nlen,NDr_idx)
      !---------------------------------------------------------------------------------
      Real, allocatable :: Events(:,:,:,:)                           ! (Nst,Ntime,Nlen,NDr_idx): =1 if Dr_IDX < Qp
      Integer :: iidx, Nyear
      Integer, Parameter :: NDr_idx=2 ! SPI, SPEI
      Integer, Parameter :: NLen=3    ! 3, 6, 12 months
      Integer, parameter :: NDr_Categ = 3             ! Moderate, Severe, Extreme

      !Real, parameter :: Fix_Thres_D(NDr_idx) = (/-0.5, -0.5/)    ! Dry or Mild drought
      !Real, parameter :: Fix_Thres_M(NDr_idx) = (/-1.0, -1.0/)    ! Moderate
      !Real, parameter :: Fix_Thres_S(NDr_idx) = (/-1.5, -1.5/)    ! Severe
      !Real, parameter :: Fix_Thres_E(NDr_idx) = (/-2.0, -2.0/)    ! Extreme

      !Real :: Fix_Thres(NDr_idx)

      Real, Parameter :: Fix_Thres(NDr_Categ)=(/-0.5,-1.0,-1.5/)
      Character (Len=3), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"spi","spe"/)
      Character (Len=3), Parameter :: Dr_Categ(NDr_Categ)=(/"mil","mod","sev"/)
                                   !(/"PDSI","SPI12","SPE12"/)
      Character (Len=250) :: Fout, Fname, LandMaskFile, RegionFile
      Character (Len=250) :: Data_Dir, Model_Obs, Out_Dir
      !---------------------------------------------------------------------------------
      Integer :: ilon,ilat,it,imon,iper,imod,isce
      Integer :: Nt0, itime1, itime2    ! Timeseries 
      !---------------------------------------------------------------------------------

      ! Drought frequency (DF): the number of events (spells)
      Real, allocatable :: DF(:,:,:)       ! (Nst,NLen,ndr_idx)
      Real, allocatable :: DFM(:,:,:,:)    ! (Nst,Nmon,Nlen,ndr_idx)
      ! Total drought duration (TDD): the sum of the durations
      Real, allocatable :: TDD(:,:,:)      ! (Nst,NLen,ndr_idx) 
      ! Total drought severity (TDS): the sum of the severities of the drought events
      Real, allocatable :: TDS(:,:,:)      ! (Nst,NLen,ndr_idx)
      ! Spatial Extent (1962-2005):
      Real, allocatable :: Geo(:,:,:)          ! (Ntime,NLen,ndr_idx)
      Real, allocatable :: GeoReg(:,:,:,:)     ! (Ntime,NLen,Nreg,ndr_idx)
      Real, allocatable :: MeanReg(:,:,:,:)    ! (Ntime,NLen,Nreg,ndr_idx)
      !---------------------------------------------------------
      Real, allocatable :: Freq_idx(:,:,:,:)     !(Nidx_bin,Nreg,Nlen,NDr_idx) idx frequency
      Real, allocatable :: Senslp(:,:,:), prob(:,:,:)        !(Nst,NLen,ndr_idx)
      !---------------------------------------------------------
      Real, allocatable :: ClimIDX(:,:)        ! (Nidx,Ntime)
      Real, allocatable :: Cor_dr_clim(:,:,:,:), pv_dr_clim(:,:,:,:) ! (Nst,Nidx,Nlen,ndr_idx)

      Real, allocatable :: eof_cli(:,:)        ! (Ntime,Nmode)
      Real, allocatable :: Cor_dr_eof_clim(:,:,:,:), pv_dr_eof_clim(:,:,:,:) ! (Nst,Nmode,Nlen,ndr_idx)
      Integer :: Nmode

      !----------------------------------------------------------------------------------------------
      ! Ket qua tinh cua chuong trinh nay (Mar 2025):

      ! DF (Nst,NLen,ndr_idx) : Trung binh nam so dot
      ! TDD(Nst,NLen,ndr_idx) : Trung binh nam so thang
      ! TDS(Nst,NLen,ndr_idx) : Trung binh nam do khac nghiet
      ! Geo(Ntime,Nlen,ndr_idx): Chuoi thoi gian ty le dien tich bi han 
      ! DFM(Nst,Nmon,Nlen,ndr_idx): Tan suat xuat hien han trong tung thang

      !----------------------------------------------------------------------------------------------
      Out_Dir = "/home/tanpv/Progs/Dr_4Books/Res/"
      !----------------------------------------------------------------------------------------------
      Fname=trim(Out_Dir)//"spi_spei_stn_1980-2020.nc"
      !Fname="/home/tanpv/Progs/Dr_VN_2020/VN_PP/Results/spi_spei_stn_1980-2020_48timescales.nc"
     
      Call Read_Drough_IDX (Fname)   ! Return Dr_IDX(nlon,nlat,Ntime,NDr_idx), xlon, xlat
      Nyear = Yr2-Yr1                ! Bo nam dau tien
      !----------------------------------------------------------------------------------------------
!      Fix_Thres(:)=Fix_Thres_D(:)    ! Drought level (D,M,S,E)
      !----------------------------------------------------------------------------------------------
!      Call Calc_Per_Events           ! if Fix_Ths=.false. then used pd3=5% --> Events(ilon,ilat,it,iidx)

      Fname="/home/tanpv/Data/SEACLID/ClimIDX_2022/Clim_IDX_2022.nc" !"/home/tanpv/Data/ClimIDX/Clim_IDX_2020.nc"
      Call Read_ClimIDX (Fname)
      Fname="/home/tanpv/Progs/Dr_4Books/Res/eof_clim.nc"
      Call Read_eof_cli (Fname)
      Call Calc_Trend
      Call Calc_corr_Dr_Clim

!      Call Calc_chars
!print*,"Calc_chars"
!      Call Calc_Trend
!print*,"Calc_Trend"
!      Call Calc_corr_Dr_Clim
!print*,"Writing..."
      ! Write out:
      ! Moderate

      Do ilev=1,NDr_Categ
         Call Calc_Per_Events (Fix_Thres(ilev))           ! if Fix_Ths=.false. then used pd3=5% --> Events(ilon,ilat,it,iidx)
         Call Calc_chars
         Do iidx=1,NDr_idx

            Fname=trim(Out_Dir)//"DF_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)      ! Moderate
            Call Write_GrADS_Stn (DF,Fname,iidx)
            Fname=trim(Out_Dir)//"TDD_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)       ! Moderate
            Call Write_GrADS_Stn (TDD,Fname,iidx)
            Fname=trim(Out_Dir)//"TDS_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)       ! Moderate
            Call Write_GrADS_Stn (TDS,Fname,iidx)
            Fname=trim(Out_Dir)//"Senslp_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)      ! Moderate
            Call Write_GrADS_Stn (Senslp,Fname,iidx)
            Fname=trim(Out_Dir)//"Prob_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)      ! Moderate
            Call Write_GrADS_Stn (prob,Fname,iidx)
            Fname=trim(Out_Dir)//"DFM_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)
            Call Write_GrADS_Stn_Multivar (DFM,Fname,iidx)
            Fname=trim(Out_Dir)//"Cor_dr_clim_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)
            Call Write_GrADS_Stn_Multivar_1 (Cor_dr_clim,Fname,iidx)
            Fname=trim(Out_Dir)//"Pval_dr_clim_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)
            Call Write_GrADS_Stn_Multivar_1 (pv_dr_clim,Fname,iidx)
            Fname=trim(Out_Dir)//"Cor_dr_eof_clim_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)
            Call Write_GrADS_Stn_Multivar_2 (Cor_dr_clim,Fname,iidx)
            Fname=trim(Out_Dir)//"Pval_dr_eof_clim_"//IDX_nm(iidx)//"_"//Dr_Categ(ilev)
            Call Write_GrADS_Stn_Multivar_2 (pv_dr_clim,Fname,iidx)
         Enddo
         Fname=trim(Out_Dir)//"geo_"//Dr_Categ(ilev)//".nc"
         Call Write_Geo(Fname)
         Fname=trim(Out_Dir)//"Freq_"//Dr_Categ(ilev)//".nc"
         Call Write_Freq_IDX (Fname)
      Enddo
   CONTAINS
    !---------------------------------------------------------------------
      Subroutine Calc_corr_Dr_Clim
         ! Cor_dr_clim(:,:,:), pv_dr_clim(:,:,:)
         ! Dr_IDX(ist,1:Ntime,ilen,1), ClimIDX(:,:)        ! (Nidx,Ntime)
         Integer :: ist, ilen, it, iidx, idr
         Real :: Cor, pval
         Real :: X(Ntime), Y(Ntime)

         Allocate (Cor_dr_clim(Nst,Nidx,Nlen,ndr_idx),pv_dr_clim(Nst,Nidx,Nlen,ndr_idx))
         Allocate (Cor_dr_eof_clim(Nst,Nmode,Nlen,ndr_idx), pv_dr_eof_clim(Nst,Nmode,Nlen,ndr_idx))
         Do idr=1,ndr_idx
            Do ilen=1,Nlen
               Do ist=1,Nst
                  X(1:Ntime) = Dr_IDX(ist,1:Ntime,ilen,idr)
                  Do iidx=1,Nidx                            ! ClimIDX
                     Y(1:Ntime) = ClimIDX(iidx,1:Ntime)
                     Call Correlation (X,Y,Ntime,Cor,pval, Miss)
                     Cor_dr_clim(ist,iidx,ilen,idr) = Cor
                     pv_dr_clim(ist,iidx,ilen,idr) = pval
                  Enddo
                  Do iidx=1,Nmode                           ! eof_cli
                     Y(1:Ntime) = eof_cli(1:Ntime,iidx)
                     Call Correlation (X,Y,Ntime,Cor,pval, Miss)
                     Cor_dr_eof_clim(ist,iidx,ilen,idr) = Cor
                     pv_dr_eof_clim(ist,iidx,ilen,idr) = pval
                  Enddo
               Enddo
            Enddo
         Enddo
      End Subroutine Calc_corr_Dr_Clim
      !--------------------------------------------------------------------------------------
      Subroutine Read_eof_cli (Fname)
       Character (Len=*) :: Fname
       Real, Allocatable :: xmod(:), xt(:)
       Integer :: Nt

       Call rd_ncdf_2d (Fname,xt,xmod,Nt,Nmode,eof_cli,"eof")
      End Subroutine Read_eof_cli
      !--------------------------------------------------------------------------------------
      Subroutine Read_ClimIDX (Fname) !,varn)
       Character (Len=*) :: Fname     !,varn
       Real, Allocatable :: Year(:),Month(:), tmp(:,:)
       Integer :: Yr11, Nyr, Nmn, iidx, iyr, i,j, k

       Allocate (ClimIDX(Nidx,Ntime))
       Yr11=Yr1+1
       Do iidx=1,Nidx
          Call rd_ncdf_2d (Fname,Year,Month,Nyr,Nmn,tmp, trim(Clim_nm(iidx)))
          j=0
          Do iyr=Yr11,Yr2
             Do i=1,Nyr
                if (Year(i)==iyr) then
                   Do k=1,Nmon
                      j=j+1; ClimIDX(iidx,j)=tmp(i,k)
                   Enddo
                endif
             Enddo
          Enddo
          Deallocate (Year,Month,tmp)
       Enddo
      End Subroutine Read_ClimIDX
   !---------------------------------------------------------------------------------
   Subroutine Calc_Trend
      Real, Allocatable :: X(:)
      Integer :: N1, ist, it, iidx,ilen, N
      Real :: Sen,pro

      Allocate (Senslp(Nst,NLen,ndr_idx),prob(Nst,NLen,ndr_idx)) !Senslp(:,:,:), prob((:,:,:)        !(Nst,NLen,ndr_idx)
      Allocate (X(Ntime))
      Senslp=Miss; prob=Miss; N1=30
      Do iidx=1,ndr_idx
         Do ilen=1,Nlen
            Do ist=1,Nst
               if (Reg(ist)>0) then
                  X(1:Ntime)=Dr_IDX(ist,1:Ntime,ilen,iidx)
                  Call MK_Sen(X,Ntime,N1,Miss,Sen,pro,N)
                  Senslp(ist,ilen,iidx)=Sen; prob(ist,ilen,iidx)=pro
               endif
            Enddo
         Enddo
      Enddo
   End Subroutine Calc_Trend
   !---------------------------------------------------------------------------------
   Subroutine Calc_chars
      Integer :: ilon,ilat,it,iidx
      Real :: X(Ntime), M1(Ntime), IMAX1, Spells(Ntime), N_Spell_1, I_Start(Ntime), I_End(Ntime)
      Integer :: N, it1,it2, i, k, N_Spell, Imax, ista, iend,imon,imon1
      Real :: tmp, ntmp,tmp1, ntmp1, Smean, Smax, In_mean, In_max, Amean, Se0, Sem, Sex, Inm, Inx

      if (Allocated (DF))  Deallocate (DF);  if (Allocated (TDD)) Deallocate (TDD)
      if (Allocated (TDS)) Deallocate (TDS); if (Allocated (Geo)) Deallocate (Geo)
      if (Allocated (GeoReg))   Deallocate (GeoReg);   if (Allocated (MeanReg)) Deallocate (MeanReg)
      if (Allocated (Freq_idx)) Deallocate (Freq_idx); if (Allocated (DFM)) Deallocate (DFM)

      Allocate (DF (Nst,NLen,ndr_idx))
      Allocate (TDD(Nst,NLen,ndr_idx))
      Allocate (TDS(Nst,NLen,ndr_idx))
      Allocate (Geo(Ntime,Nlen,ndr_idx))
      Allocate (GeoReg(Ntime,Nlen,Nreg,ndr_idx))
      Allocate (MeanReg(Ntime,Nlen,Nreg,ndr_idx))
      Allocate (Freq_idx(Nidx_bin,Nreg,Nlen,NDr_idx))
      Allocate (DFM(Nst,Nmon,Nlen,ndr_idx))

      print*,"Calculating_chars"
      DF=Miss; TDD=Miss; TDS=Miss; Geo=Miss; GeoReg=Miss
      M1=Miss; Spells=Miss; I_Start=Miss;I_End=Miss
      MeanReg=Miss; DFM=Miss

      ! DF (Nst,NLen,ndr_idx) : Trung binh nam so dot
      ! TDD(Nst,NLen,ndr_idx) : Trung binh nam so thang
      ! TDS(Nst,NLen,ndr_idx) : Trung binh nam do khac nghiet
      ! Geo(Ntime,Nlen,ndr_idx): Chuoi thoi gian ty le dien tich bi han 
      ! DFM(Nst,Nmon,Nlen,ndr_idx): Tan suat xuat hien han trong tung thang

      Do iidx=1,NDr_idx
         Do ilen=1, Nlen
            Do ist=1, Nst
               if (Reg(ist) > 0) then
                  X(:) = Events(ist,:,ilen,iidx)  ! Ntime values
                  Call NSpell_and_Leng (X, Ntime, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss)
                  N_Spell=Int(N_Spell_1)
                  DF (ist,ilen,iidx) = N_Spell_1/Real(Nyear)  ! per year - Annual average
                  tmp=0
                  Do i=1,N_Spell
                     tmp=tmp+Spells(i)    ! Tong do dai tat ca cac dot
                  Enddo
                  TDD(ist,ilen,iidx) = tmp/Real(Nyear) ! per year - Trung binh nam So thang co han  
                  tmp=0
                  Do it=1,Ntime
                     if (X(it)==1.AND.Dr_IDX(ist,it,ilen,iidx)/=Miss) tmp = tmp + Dr_IDX(ist,it,ilen,iidx) 
                  Enddo
                  TDS(ist,ilen,iidx) = tmp/Real(Nyear)  ! per year - Trung binh nam severity
                  Do imon=1,Nmon       ! Monthly drought(DFM)
                     tmp=0; ntmp=0
                     Do it=1,Ntime
                        imon1=mod(it,Nmon); if (imon1==0) imon1=Nmon
                        if (imon1==imon) then
                           ntmp=ntmp+1
                           if (Events(ist,it,ilen,iidx)==1) tmp=tmp+1
                        endif
                     Enddo
                     DFM(ist,imon,ilen,iidx)=tmp/ntmp*100.0
                  Enddo
               endif
            Enddo
         Enddo
         ! Spatial Extent:
         Do ilen=1,Nlen
            Do it=1,Ntime
               tmp=0; ntmp=0
               Do ist=1, Nst
                  if (Reg(ist) >0 ) then
                     ntmp=ntmp+1                                  ! Number of landpoints
                     if (Events(ist,it,ilen,iidx)==1) tmp=tmp+1  ! Number of drought points
                  endif
               Enddo
               if (ntmp/=0) Geo(it,ilen,iidx) = tmp/ntmp*100              ! Percentage
            Enddo
         Enddo
         ! GeoReg:
         Do ilen=1,Nlen
            Do ireg=1,Nreg
               Do it=1, Ntime
                  tmp=0; ntmp=0
                  Do ist=1, Nst
                     if (Reg(ist)==ireg) then
                        ntmp=ntmp+1
                        if (Events(ist,it,ilen,iidx)==1) tmp=tmp+1
                     endif
                  Enddo
                  if (ntmp/=0) GeoReg(it,ilen,ireg,iidx) = tmp/ntmp*100
               Enddo
            Enddo
         Enddo
         ! MeanReg:
         Do ilen=1,Nlen
            Do ireg=1,Nreg
               Do it=1, Ntime
                  tmp=0; ntmp=0
                  Do ist=1, Nst
                     if (Reg(ist)==ireg) then
                        ntmp=ntmp+1
                        tmp=tmp+Dr_IDX(ist,it,ilen,iidx) 
                     endif
                  Enddo
                  if (ntmp/=0) MeanReg(it,ilen,ireg,iidx) = tmp/ntmp
               Enddo
            Enddo
         Enddo
      Enddo 
      ! IDX frequency:
      ! Freq_idx(Nidx_bin,Nreg,NDr_idx)) 
      Do iidx=1,NDr_idx
         Do ireg=1,Nreg
            Do ilen=1,Nlen
               Do ibin=1,Nidx_bin
                  tmp=0; ntmp=0        ! Ca giai doan
                  Do ist=1, Nst
                     if (Reg(ist)==ireg) then
                        Do it=1,Ntime
                           if (Dr_IDX(ist,it,ilen,iidx)/=Miss) then
                              tmp=tmp+1
                              if (Dr_IDX(ist,it,ilen,iidx)>=Idx_bin(ibin-1).AND.Dr_IDX(ist,it,ilen,iidx)<Idx_bin(ibin)) ntmp=ntmp+1
                           endif
                        Enddo
                     endif
                  Enddo
                  Freq_idx(ibin,ireg,ilen,iidx)=ntmp/tmp
               Enddo
            Enddo
         Enddo
      Enddo
   End Subroutine Calc_chars
   !---------------------------------------------------------------------------------
   Subroutine Write_Freq_IDX (Fname)
      Character (Len=*) :: Fname
      Character (Len=50)  :: var_n, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=250) :: var_ln
      Real :: x_idx(NDr_idx), xreg(Nreg), xlen(Nlen)
      Integer :: i,j,k,ilen

      ! Freq_idx(:,:,:)     !(Nidx_bin,Nreg,NDr_idx)
      Do i=1, NDr_idx; x_idx(i)=i; Enddo
      Do i=1,Nreg; xreg(i)=i; Enddo
      Do ilen=1,Nlen; xlen(ilen)=ilen; Enddo

      lon_nm="bin"; lat_nm="reg"; z_nm="len"; t_nm="idx"; var_unit="-"
      t_unit="years since 1981-01-01 "
      var_ln="Idx frequency of spi and spei for different regions"
      Call wr_ncdf_4d (Fname,Xidx_bin,xreg,xlen,x_idx,Nidx_bin,Nreg,Nlen,NDr_idx,Freq_idx,"freq", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)

      !Call wr_ncdf_3d (Fname,Xidx_bin,xreg,x_idx,Nidx_bin,Nreg,NDr_idx,Freq_idx,"freq",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss)
      
      ! wr_ncdf_3d (filenm,lon,lat,t,nlon,nlat,nt,var,var_nm,t_unit, var_long_name, var_unit, lon_nm, lat_nm, t_nm, flag)
   End Subroutine Write_Freq_IDX
   !---------------------------------------------------------------------------------
   Subroutine Write_Geo(Fname)
      Character (Len=*) :: Fname
      Character (Len=50)  :: var_n, var_unit, t_unit, lon_nm, lat_nm, z_nm, t_nm
      Character (Len=250) :: var_ln
      Real :: x_idx(NDr_idx), xmon(Ntime), xreg(Nreg), xlen(Nlen), xreg1(Nreg+1)
      Real, Allocatable, Dimension (:,:,:,:) :: X4d, Y4d
      Integer :: ireg, iidx, iyr

      lon_nm="time"; lat_nm="len"; z_nm="idx"; t_nm="reg";
      t_unit="years since 1981-01-01 "
      Do iidx=1,NDr_idx; x_idx(iidx)=iidx; Enddo
      Do it=1,Ntime; xmon(it)=it; Enddo
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Do ilen=1,Nlen; xlen(ilen)=ilen; Enddo
      
      ! Geographic Extent:
      ! Regions: 
      Allocate (X4d(Ntime,NLen,ndr_idx,Nreg),Y4d(Ntime,NLen,ndr_idx,Nreg))
      Do it=1,Ntime
         Do ilen=1,Nlen
            Do ireg=1,Nreg
               Do iidx=1,ndr_idx
                  X4d(it,ilen,iidx,ireg)=GeoReg (it,ilen,ireg,iidx)
                  Y4d(it,ilen,iidx,ireg)=MeanReg(it,ilen,ireg,iidx)
               Enddo
            Enddo
         Enddo
      Enddo
      var_ln="Spatial drought extent for different timescales and regions using SPI and SPEI"
      var_unit="%"
      Call wr_ncdf_4d (Fname,xmon,xlen,x_idx,xreg,Ntime,Nlen,NDr_idx,Nreg,X4d,"georeg", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, z_nm, t_nm, Miss)
      ! MeanReg:
      var_ln="Mean values of drought index"
      var_unit="-"
      Call wr_ncdf_4d_apd (Fname,Y4d,"mean",Ntime,Nlen,NDr_idx,Nreg, &
                           var_ln,var_unit, lon_nm, lat_nm, z_nm, t_nm,Miss)
      ! VN: Geo(:,:,:)          ! (Ntime,NLen,ndr_idx)
      var_ln="Spatial drought extent for different timescales over entire country using SPI and SPEI"
      Call wr_ncdf_3d_apd (Fname, Geo,"geovn",Ntime,Nlen,ndr_idx, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      Deallocate (X4d,Y4d)
   End Subroutine Write_Geo
   !---------------------------------------------------------------------------------
   Subroutine Calc_Per_Events (Thres)
      ! Return: Qp and Events
      Integer :: ilon, ilat, iidx, it
      Real :: p,Thres
      Real :: X(Ntime)
         
      if (Allocated (Events)) Deallocate (Events)
      Allocate (Events(Nst,Ntime,Nlen,NDr_idx))
      Events = Miss;  ! Need to be check !!!

      Do iidx=1,NDr_idx
         Do ilen=1, Nlen
            Do ist=1, Nst
               if (Reg(ist) > 0) then
                  Do it=1, Ntime
                     if (Dr_IDX(ist,it,ilen,iidx) /=Miss .AND. Dr_IDX(ist,it,ilen,iidx) <= Thres) then ! Fix_Thres(iidx)) then 
                        Events(ist,it,ilen,iidx) = 1
                     else
                        Events(ist,it,ilen,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, Nle
         Character (Len=*) :: Fname
         Real :: tmp, xmean(12), ic(12)
         Real, Allocatable :: Stn(:), time(:), Xlen(:), X1d(:), Y1d(:)

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

         Call rd_ncdf_3d (trim(Fname), Stn, time, xlen, Nst, Nt, nle, SPI, "spi")
         Deallocate (Stn, time, xlen)
         Call rd_ncdf_3d (trim(Fname), Stn, time, xlen, Nst, Nt, nle, SPEI, "spei") 
         Deallocate (Stn, time, xlen)

         Call rd_ncdf_1d (Fname, Stn, Nst, xlon, "lon")
         Deallocate (Stn)
         Call rd_ncdf_1d (Fname, Stn, Nst, xlat, "lat")
         Deallocate (Stn)
         Call rd_ncdf_1d (Fname, Stn, Nst, Reg, "reg")
         Deallocate (Stn)
         Call rd_ncdf_1d_Char (Fname, Stn, Nst, Stname, "stname")
         Deallocate (Stn)

         Where (IsNaN(SPEI)); SPEI=Miss; Endwhere
         Where (IsNaN(SPI)); SPI=Miss; Endwhere

         Ntime = Nt-Nmon     ! Skip first year due to SPI and SPEI were missed
         it1 = Nmon+1; it2 = Nt

         Allocate (Dr_IDX(Nst,Ntime,Nlen,NDr_idx), X1d(Ntime), Y1d(Ntime))

         Dr_IDX(:,:,:,:) = Miss

         Do ist=1, Nst
            Do ilen=1,Nlen
               X1d(1:Ntime)=SPI (ist,it1:it2,ilen)
               Y1d(1:Ntime)=SPEI(ist,it1:it2,ilen)
               if (.not.Empty_Ar1(X1d, Ntime).AND. &
                   .not.Empty_Ar1(Y1d, Ntime) ) then 
                  Dr_IDX(ist,1:Ntime,ilen,1)=X1d(1:Ntime)  ! SPI
                  Dr_IDX(ist,1:Ntime,ilen,2)=Y1d(1:Ntime)  ! SPEI
               else
                  Reg(ist)=Miss; exit
               endif
            Enddo
         Enddo
         Deallocate (X1d,Y1d)
         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/Progs/Dr_4Books/Sub_Write_to_GrADS.f90"
   Include "/home/tanpv/Mylibs/Calc_Correlation_Coefficient.f90"
   !---------------------------------------------------------------------------------
   END Program
   Include "/home/tanpv/Mylibs/Student_Prob.f90"