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

      ! 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(:,:,:,:)                           ! (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 :: GeoReg(:,:,:)                             ! (Ntime,NDr_idx,Nreg) - 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

      ! New define for Severity and Intensity:
      ! Severity = Sum of DrIDX-Threshold; Intensity = Severity/Duration:
      !Real, allocatable :: SEVm(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Mean Severity
      !Real, allocatable :: SEVx(:,:,:)                              ! ((Nlon,Nlat,NDr_idx) - Max Severity
      !Real, allocatable :: INTEm(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Mean Intensity
      !Real, allocatable :: INTEx(:,:,:)                            ! ((Nlon,Nlat,NDr_idx) - Max Intensity


      !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 :: Geo_PDF(:,:,:)                            ! (NDr_idx,Nreg,NGeo_bin) PDF
      !Real, allocatable :: GeoRegMon(:,:,:)                          ! (Ntime,NDr_idx,Nreg+1) - Geo Ext in each mon

      Integer :: Nst, iidx, Nyear
      Integer, Parameter :: NDr_idx=3 ! PDSI, 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.,.true./)     ! Use Fix_Thres if Fix_Ths=.true.

      Character (Len=5), Parameter :: IDX_nm(NDr_idx)= &
                                      (/"pdsi","spi12","spe12"/)
                                   !(/"PDSI","SPI12","SPE12"/)
      Integer, Parameter :: NPeriod=3, NSce=2, Nmodels=8 ! 8=6-models+Ens+Obs
      Character (Len=250) :: Fout, Fname, LandMaskFile, RegionFile
      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, 1962-2005 (5 periods)
      ! Drought frequency (DF): the number of events (spells)
      Integer, parameter :: Ndec=5
      Real, allocatable :: DF(:,:,:,:)       ! (nlon,nlat,ndec,ndr_idx)
      Real, allocatable :: DFY(:,:,:,:)      ! (nlon,nlat,nyear,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

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

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

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

      Nyear = Yr2-Yr1           ! Bo nam dau tien
      Nyrsec = Nyear/(Ndec-1)   ! Moi thap ky (11 nam)
      
      Fname="/home/tanpv/Progs/Dr_SEA/Res_APH/BL_APH.nc"
      Call Read_Drough_IDX (Fname)    ! Return Dr_IDX(nlon,nlat,Ntime,NDr_idx), xlon, xlat

      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)   ! =528/4 = 132

         Call Calc_chars
         Call Calc_Trend_DFY_TDSY_GeoReg

         ! Write out:
         Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Chars_"//Sevnm(isev)//".nc"      ! Moderate
         Fout ="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_"//Sevnm(isev)//".nc"
         RegionFile="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_reg_"//Sevnm(isev)//".nc"

      !Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Chars_Sev.nc"   ! Severe
      !Fout ="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_Sev.nc"
      !RegionFile="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_reg_Sev.nc"

      !Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Chars_Ext.nc"  ! Extreme
      !Fout ="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_Ext.nc"
      !RegionFile="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_geo_reg_Ext.nc"

         Call Write_Char_to_NC (Fname, Fout,RegionFile)
         Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Freq_Idx_"//Sevnm(isev)//".nc"
         Call Write_Freq_IDX (Fname)
         Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Areal_mean_"//Sevnm(isev)//".nc"
         Call Write_Areal_mean_of_drought (Fname)
         Fname="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Trend_DFY_TDSY_"//Sevnm(isev)//".nc"
         RegionFile="/home/tanpv/Progs/Dr_SEA/Calc_Chars/Out_APH/APH_Trend_GeoReg_"//Sevnm(isev)//".nc"
         Call Write_Trend_DFY_TDSY_to_NC(Fname,RegionFile)

         Deallocate (Qp,Events)
         Deallocate (DF, TDD, TDS, DFY, TDSY, 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 :: N, it1,it2, i, k, N_Spell, Imax, ista, iend
      Real :: tmp, ntmp,tmp1, ntmp1, Smean, Smax, In_mean, In_max, Amean, Se0, Sem, Sex, Inm, Inx

      Allocate (DF (nlon,nlat,ndec,ndr_idx))
      Allocate (TDD(nlon,nlat,ndec,ndr_idx))
      Allocate (TDS(nlon,nlat,ndec,ndr_idx))

      Allocate (DFY (nlon,nlat,Nyear,ndr_idx))
      !Allocate (TDDY(nlon,nlat,Nyear,ndr_idx))
      Allocate (TDSY(nlon,nlat,Nyear,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

      Do iidx=1,NDr_idx
         Do ilon=1, Nlon
            Do ilat=1, Nlat
               if (Mask(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)
                  N_Spell=Int(N_Spell_1)
                  DF (ilon,ilat,ndec,iidx) = N_Spell_1/Nyear*10.  ! per decade               ! Annual average
                  TDD(ilon,ilat,ndec,iidx) = Sum(Spells(1:N_Spell))/Nyear*10.  ! per decade 
                  tmp=0
                  Do it=1,Ntime
                     if (X(it)==1.AND.Dr_IDX(ilon,ilat,it,iidx)/=Miss) tmp = tmp + Dr_IDX(ilon,ilat,it,iidx)
                  Enddo
                  TDS(ilon,ilat,ndec,iidx) = tmp/Nyear*10.  ! per decade

                  Do k=1,ndec-1
                     it1=(k-1)*nsec+1; it2=it1+nsec-1
                     X1(1:Nsec) = Events(ilon,ilat,it1:it2,iidx)
                     Call NSpell_and_Leng (X1, Nsec, M1, Imax1, Spells, N_Spell_1, I_Start, I_End, Miss)
                     N_Spell=Int(N_Spell_1)
                     DF (ilon,ilat,k,iidx) = N_Spell_1/Nyrsec*10.  ! per decade               ! Annual average
                     TDD(ilon,ilat,k,iidx) = Sum(Spells(1:N_Spell))/Nyrsec*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(ilon,ilat,it,iidx)/=Miss) tmp = tmp + Dr_IDX(ilon,ilat,it,iidx)
                     Enddo
                     TDS(ilon,ilat,k,iidx) = tmp/Nyrsec*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(ilon,ilat,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(ilon,ilat,it,iidx)/=Miss) then
                           tmp=tmp+1; tmp1=tmp1+Dr_IDX(ilon,ilat,it,iidx)      ! So thang han va tong Dr_IDX<Thres
                        endif
                     Enddo
                     DFY (ilon,ilat,iyr,iidx) = tmp; TDSY(ilon,ilat,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
         ! Spatial Extent:
         Do it=1,Ntime
            tmp=0; ntmp=0
            Do ilon=1, Nlon
               Do ilat=1, Nlat
                  if (Mask(ilon,ilat) /= Miss) then
                     ntmp=ntmp+1                                  ! Number of landpoints
                     if (Events(ilon,ilat,it,iidx)==1) tmp=tmp+1  ! Number of drought points
                  endif
               Enddo
            Enddo
            if (ntmp/=0) Geo(it,iidx) = tmp/ntmp*100              ! Percentage
            ! Reg:
            Do ireg=1,Nreg
               tmp=0; ntmp=0
               Do ilon=1, Nlon
                   Do ilat=1, Nlat
                      if (Mask(ilon,ilat) /= Miss.AND.Regs(ilon,ilat)==ireg) then
                         ntmp=ntmp+1
                         if (Events(ilon,ilat,it,iidx)==1) tmp=tmp+1
                      endif
                   Enddo
               Enddo
               if (ntmp/=0) GeoReg(it,ireg,iidx) = tmp/ntmp*100
            Enddo
         Enddo
      Enddo 
      ! IDX frequency:
      ! Freq_idx(:,:,:)     !(Nidx_bin,Ndec,NDr_idx) idx frequency
      Do iidx=1,NDr_idx
         Do ibin=1,Nidx_bin
            tmp=0; ntmp=0        ! Ca giai doan
            Do ilon=1, Nlon
               Do ilat=1, Nlat
                  if (Mask(ilon,ilat) /= Miss) then
                     Do it=1,Ntime
                        if (Dr_IDX(ilon,ilat,it,iidx)/=Miss) then
                           tmp=tmp+1
                           if (Dr_IDX(ilon,ilat,it,iidx)>=Idx_bin(ibin-1).AND.Dr_IDX(ilon,ilat,it,iidx)<Idx_bin(ibin)) ntmp=ntmp+1
                        endif
                     Enddo
                  endif
               Enddo
            Enddo
            Freq_idx(ibin,Ndec,iidx)=ntmp/tmp
            Do k=1,ndec-1       ! Tung thap ky
               tmp=0; ntmp=0
               it1=(k-1)*nsec+1; it2=it1+nsec-1
               Do ilon=1, Nlon
                  Do ilat=1, Nlat
                     if (Mask(ilon,ilat) /= Miss) then
                        Do it=it1,it2   ! So thang trong thap ky
                           if (Dr_IDX(ilon,ilat,it,iidx)/=Miss) then
                              tmp=tmp+1
                              if (Dr_IDX(ilon,ilat,it,iidx)>=Idx_bin(ibin-1).AND.Dr_IDX(ilon,ilat,it,iidx)<Idx_bin(ibin)) ntmp=ntmp+1
                           endif 
                        Enddo
                     endif
                  Enddo
               Enddo
               Freq_idx(ibin,k,iidx)=ntmp/tmp
            Enddo
         Enddo
      Enddo
      ! Areal mean:
      Do iidx=1,NDr_idx
         Do it=1,Ntime
            tmp=0; ntmp=0
            Do ilon=1, Nlon
               Do ilat=1, Nlat
                  if (Mask(ilon,ilat) /= Miss) then
                     if (Dr_IDX(ilon,ilat,it,iidx)/=Miss) then
                        ntmp=ntmp+1; tmp=tmp+Dr_IDX(ilon,ilat,it,iidx)
                     endif
                  endif
               Enddo
            Enddo
            Areal_mean(it,iidx)=tmp/ntmp
         Enddo
      Enddo

   End Subroutine Calc_chars
   !---------------------------------------------------------------------------------
   Subroutine Calc_Trend_DFY_TDSY_GeoReg
        Integer :: i,j,NumYear,N1,ilon,ilat
        Real :: SenSlope, prob, X1d(Nyear), Y1d(Ntime), Ptb

        NumYear=20    ! Minimum sample
        print*,"Calculating Trend drought..."

        Allocate (Sen_DFY(nlon,nlat,ndr_idx),  prob_DFY(nlon,nlat,ndr_idx))
        Allocate (Sen_TDSY(nlon,nlat,ndr_idx), prob_TDSY(nlon,nlat,ndr_idx))
        Allocate (Sen_GeoReg(Nreg,ndr_idx), prob_GeoReg(Nreg,ndr_idx))

        Sen_DFY=Miss; Sen_TDSY=Miss; Sen_GeoReg=Miss
        prob_DFY=Miss; prob_TDSY=Miss; prob_GeoReg=Miss

        Do iidx=1,Ndr_idx
           Do ilon=1,Nlon
              Do ilat=1,Nlat
                 if (Mask(ilon,ilat) /= Miss) then
                    X1d(1:Nyear) = DFY(ilon,ilat,1:Nyear,iidx)
                    Call MK_Sen(X1d,Nyear,10,Miss,SenSlope,prob,N1)
                    Sen_DFY(ilon,ilat,iidx)=SenSlope; prob_DFY(ilon,ilat,iidx)=prob
                    X1d(1:Nyear) = TDSY(ilon,ilat,1:Nyear,iidx)
                    Call MK_Sen(X1d,Nyear,10,Miss,SenSlope,prob,N1)
                    Sen_TDSY(ilon,ilat,iidx)=SenSlope; prob_TDSY(ilon,ilat,iidx)=prob
                 endif
              Enddo
           Enddo
        Enddo
        !
        Do iidx=1,Ndr_idx
           Do ireg=1,Nreg
              Y1d(1:Ntime) = GeoReg(1:Ntime,ireg,iidx)
              Call MK_Sen(Y1d,Ntime,20,Miss,SenSlope,prob,N1)
              Sen_GeoReg(ireg,iidx)=SenSlope; prob_GeoReg(ireg,iidx)=prob
           Enddo
        Enddo
   End Subroutine Calc_Trend_DFY_TDSY_GeoReg
   !---------------------------------------------------------------------------------
   Subroutine Write_Trend_DFY_TDSY_to_NC(Fname,Fn_geo)
      Character (Len=*) :: Fname, Fn_geo
      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)
      Real, Allocatable, Dimension (:,:,:) :: Xtmp
      Integer :: ireg, iidx, iyr

      lon_nm="lon"; lat_nm="lat"; t_nm="idx";
      t_unit="years since 1961-01-01 "
      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3
      ! Trend of DFY_TDSY 
      var_ln="Trend of yearly number of drought months (mon/year)"
      var_unit="mon/year"
      Call wr_ncdf_3d (Fname,xlon,xlat,x_idx,nlon,nlat,Ndr_idx,Sen_DFY,"dfy", &
             t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Probability of Trend of yearly number of drought months (-)"; var_unit="-"
      Call wr_ncdf_3d_apd (Fname,prob_DFY,"p_dfy",nlon,nlat,Ndr_idx, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Trend of yearly drought severity (1/year)"
      var_unit="1/year"
      Call wr_ncdf_3d_apd (Fname,Sen_TDSY,"tdsy",nlon,nlat,Ndr_idx, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss)
      var_ln="Probability of Trend of yearly drought severity"
      var_unit="-"
      Call wr_ncdf_3d_apd (Fname,prob_TDSY,"p_tdsy",nlon,nlat,Ndr_idx, var_ln,var_unit, lon_nm, lat_nm, t_nm, Miss)

      ! Trend of GeoReg:
      Do ireg=1, Nreg; xreg(ireg)=ireg; Enddo
      var_ln="Trend of geographic extend for each region (%/month)"
      var_unit="%/month"; lon_nm="reg"; lat_nm="idx";
      Call wr_ncdf_2d (Fn_geo,xreg,x_idx,Nreg,NDr_idx,Sen_GeoReg,"geo", var_ln, var_unit, lon_nm, lat_nm, Miss)
      var_ln="Probability of Trend of geographic extend for each region"
      var_unit="-"
      Call wr_ncdf_2d_apd (Fn_geo, prob_GeoReg,"p_geo",Nreg,NDr_idx,var_ln,var_unit, lon_nm, lat_nm, Miss)
   End Subroutine Write_Trend_DFY_TDSY_to_NC
   !---------------------------------------------------------------------------------
   Subroutine Read_Land_Mask
      Integer :: nx, ny
      Real, Allocatable :: lonx(:), latx(:)

      !LandMaskFile="/home/tanpv/Data/AWC/Land_Mask_VN.nc"
      !Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"all") ! all=2..8=7 regs

      Call rd_ncdf_2d (trim(LandMaskFile), xlon,xlat, nlon,nlat,Mask,"land") ! all=2..8=7 regs
      !Call rd_ncdf_2d (trim(RegionFile), lonx, latx, nx, ny, Regs, "reg")
  
      !Allocate (Mask_all(nlon,nlat))
      !Mask_all=Mask
   End Subroutine Read_Land_Mask
   !---------------------------------------------------------------------------------
   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), xdec(Ndec)
      Integer :: i,j,k

      xdec(1)=1962; xdec(2)=1973; xdec(3)=1984; xdec(4)=1995; xdec(5)=2005
      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3

      lon_nm="bin"; lat_nm="dec"; t_nm="idx"; var_unit="-"
      t_unit="years since 1961-01-01 "
      var_ln="Idx frequency: IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_3d (Fname,Xidx_bin,xdec,x_idx,Nidx_bin,Ndec,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_Char_to_NC(Fname,Fn_geo, Fn_georeg)
      Character (Len=*) :: Fname, Fn_geo, Fn_georeg
      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), xdec(Ndec), xreg1(Nreg+1)
      Real, Allocatable, Dimension (:,:,:) :: Xtmp
      Integer :: ireg, iidx, iyr

      lon_nm="lon"; lat_nm="lat"; t_nm="idx";
      t_unit="years since 1961-01-01 "
      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3
      xdec(1)=1962; xdec(2)=1973; xdec(3)=1984; xdec(4)=1995; xdec(5)=2005
      ! Drought duration:
      var_ln="Annual frequency (DF) for 62-72,73-83,84-94,95-05 and 1962-2005. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="spell"
      Call wr_ncdf_4d (Fname,xlon,xlat,xdec,x_idx,nlon,nlat,Ndec,NDr_idx,DF,"df", &
                       t_unit, var_ln, var_unit, lon_nm, lat_nm, "decade", "idx", Miss)

      var_ln="Annual drought duration (TDD) for 62-72,73-83,84-94,95-05 and 1962-2005. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="mon"
      Call wr_ncdf_4d_apd (Fname, TDD, "tdd", nlon,nlat,Ndec,NDr_idx, &
                           var_ln, var_unit, lon_nm, lat_nm, "decade", "idx", Miss)

      var_ln="Annual drought drought severity (TDS) for 62-72,73-83,84-94,95-05 and 1962-2005. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      var_unit="-"
      Call wr_ncdf_4d_apd (Fname, TDS, "tds", nlon,nlat,Ndec,NDr_idx, &
                           var_ln, var_unit, lon_nm, lat_nm, "decade", "idx", Miss)


      ! Geographic Extent: Geo(:,:) ! (Ntime,NDr_idx)
      Do iyr=1, Ntime; xmon(iyr)=iyr; Enddo
      var_ln="Timeseries (month) of Drought geographic extent (%) since Jan 1962. 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)

      var_ln="Timeseries (month) of Regional drought geographic extent (%) since Jan 1962. IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      lon_nm="mon"; lat_nm="reg"; t_nm="idx";
      Do ireg=1,Nreg; xreg(ireg)=ireg; Enddo
      Call wr_ncdf_3d (Fn_georeg,xmon,xreg,x_idx,Ntime,Nreg,NDr_idx,GeoReg,"georeg",t_unit, var_ln, var_unit, lon_nm, lat_nm, t_nm,Miss)

   End Subroutine Write_Char_to_NC
   !---------------------------------------------------------------------------------
   Subroutine Write_Areal_mean_of_drought (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)
      Real, Allocatable, Dimension (:,:,:) :: Xtmp
      Integer :: i, iidx, iyr

      x_idx(1)=1; x_idx(2)=2; x_idx(3)=3
      Do i=1,Ntime; xmon(i)=i; Enddo

      lon_nm="time"; lat_nm="idx"; var_unit="-"
      var_ln="Timeseries of areal mean of drought index: IDXs is 1:PDSI, 2:SPI12, 3:SPEI12"
      Call wr_ncdf_2d (Fname,xmon,x_idx,Ntime,NDr_idx,Areal_mean,"mean", var_ln, var_unit, lon_nm, lat_nm, Miss)

   End Subroutine Write_Areal_mean_of_drought
   !---------------------------------------------------------------------------------
   Subroutine Calc_Per_Events (p)
      ! Return: Qp and Events
      Integer :: ilon, ilat, iidx, it
      Real :: p
      Real :: X(Ntime)
         
      Allocate (Qp(Nlon,Nlat,NDr_idx),Events(Nlon,Nlat,Ntime,NDr_idx))

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

      Do iidx=1,NDr_idx
         Do 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(ilon,ilat,iidx) = Fix_Thres(iidx)
                  else
                     X(:) = Dr_IDX(ilon,ilat,:,iidx)
                     Qp(ilon,ilat,iidx) = percentile (X,Ntime,p,Miss)
                  endif
                  Do it=1, Ntime
                     if (Dr_IDX(ilon,ilat,it,iidx) /=Miss .AND. 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
               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, "spei_12")
         Deallocate (xlon, xlat, xtime)
         Call rd_ncdf_3d (trim(Fname), xlon, xlat, xtime, nlon, nlat, nt, SPI12, "spi_12")
         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

         Allocate (X1(nlon, nlat, Ntime),X2(nlon, nlat, Ntime),X3(nlon, nlat, Ntime))

         X1(:,:,1:Ntime)=PDSI  (:,:,it1:it2) !Nmon+1:Nt)
         X2(:,:,1:Ntime)=SPEI12(:,:,it1:it2) !Nmon+1:Nt)
         X3(:,:,1:Ntime)=SPI12 (:,:,it1:it2) !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

         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(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