use netcdf
Implicit None
Real, Parameter :: Miss = -99.
Integer, Parameter :: Nreg = 7, Nmon = 12, Npar = 42  ! 6 Trung binh cua cac chi so tren 7 vung KH
Integer, Parameter :: Nidx_Area = 144       ! So truong hop cac chi so tinh dien tich han (Moderate, Severe & Extreme cua PDSI, scPDSI, SPI01-12
Real, Allocatable, Dimension (:,:) :: PDSI, WPLM, ZIND, scPDSI, scWPLM, scZIND
Real, Allocatable, Dimension (:,:) :: SPI01, SPI03, SPI06, SPI12
Real, Allocatable, Dimension (:,:) :: Precip
Real, Allocatable, Dimension (:)   :: Reg
Real, Allocatable, Dimension (:)     :: Xlon, Xlat, Xtime
!Real, Allocatable, Dimension (:)     :: tmplon, tmplat, tmpTime

Integer :: ivar, ist, it, Ntime, Nst
!Integer :: Nlatmp, Nlotmp, Ntmp

Character (Len=255) :: FName
Character (Len=20)  :: VarName, Units

Logical :: status
!---------------------------------------------
Real, Allocatable, Dimension (:) :: TQ_PDSI_SPI01, TQ_PDSI_SPI03, TQ_PDSI_SPI06, TQ_PDSI_SPI12  ! Corr
Real, Allocatable, Dimension (:) :: Pr_PDSI_SPI01, Pr_PDSI_SPI03, Pr_PDSI_SPI06, Pr_PDSI_SPI12  ! Prob
Real, Allocatable, Dimension (:) :: TQ_Precip_SPI01, TQ_Precip_SPI03, TQ_Precip_SPI06, TQ_Precip_SPI12
! Area Means: 
Real, Allocatable, Dimension (:,:) :: Ave_IDX  ! (1..42,1:Ntime): PDSI, scPDSI, SPI01-12
! Geographical Extent (theo nguong):
Real, Allocatable, Dimension (:,:) :: Percent_Area
! Drought events: 0/1, depending on drought severity (Nguong)
Real, Allocatable, Dimension (:,:) :: Event_Dr
Real, Allocatable, Dimension (:,:,:) :: Events, Value_IDX
! Do dai trung binh, cuc dai cac dot han; Khoang thoi gian trung binh, cuc dai giua hai dot han:
Real, Allocatable, Dimension (:,:) :: Spell_Ave, Spell_Max, Time_Arr_Ave, Time_Arr_Max

! Trend of 6 indices of PDSI, scPDSI, SPI01, SPI03, SPI06, SPI12
! MK_Sen(Y,N,NumYear,Miss,SenSlope,prob,N1)
Real, Allocatable, Dimension (:,:) :: SenIDX, ProbIDX    ! Sen's slopes & Probability of indices

!--------------------------------------------------------------------------------
! Reading PDSI
   Fname="/home/tanpv/SEACLID_New/RES/PDSI_Res/PDSI_Stn.nc"
   Call Read_PDSI_nc (Fname)
! Reading SPI (SPI06, 12--> SPI06,12)
!   Fname="/home/tanpv/SEACLID_New/RES/SPI_Res/SPI_Stn.nc"
!   Call Read_SPI_nc (FName)
! Reading SPEI (SPEI06,12 -->SPI01,03)
   Fname="/home/tanpv/SEACLID_New/RES/SPI_Res/SPEI_Stn.nc"
   Call Read_SPEI_nc (FName)
! Reading SPI (SPI06, 12--> SPI06,12)
   Fname="/home/tanpv/SEACLID_New/RES/SPI_Res/SPI_Stn.nc"
   Call Read_SPI_nc (FName)

! Write out IDX to Text file
Call Write_IDX_To_TEXT (PDSI, Nst, Ntime, xLon, xLat, Reg, "PDSI", "The PDSI drought index", "PDSI.txt")
Call Write_IDX_To_TEXT (scPDSI, Nst, Ntime, xLon, xLat, Reg, "scPDSI", "The scPDSI drought index", "scPDSI.txt")
Call Write_IDX_To_TEXT (SPI01, Nst, Ntime, xLon, xLat, Reg, "SPEI06", "The SPEI-06 drought index", "SPEI_06.txt")
Call Write_IDX_To_TEXT (SPI03, Nst, Ntime, xLon, xLat, Reg, "SPEI12", "The SPEI-12 drought index", "SPEI_12.txt")
Call Write_IDX_To_TEXT (SPI06, Nst, Ntime, xLon, xLat, Reg, "SPI06", "The SPI-06 drought index", "SPI_06.txt")
Call Write_IDX_To_TEXT (SPI12, Nst, Ntime, xLon, xLat, Reg, "SPI12", "The SPI-12 drought index", "SPI_12.txt")

!>>>>>>>>>>>>>>> CHU Y: SPI01=SPEI06, SPI03=SPEI12
!================================================================================
   Call TQ_PDSI_SPI            ! Write out "Tq_pdsi_spi.nc"
   print*,"Call TQ_PDSI_SPI completed ..."
   Call Ave_subreg_PDSI_SPI    ! Write out "Ave_Drought_IDX.nc", "Percent_Dr_Area.nc" 
   print*,"Call Ave_subreg_PDSI_SPI completed ..."
   Call Drought_Events         ! Write out "Drought_Events.nc", "Num_Drou_Mon.nc", "Fre_Mag.nc"
   print*,"Call Drought_Events completed ..."
   Call Cal_Trend_Sen
   print*,"Call Cal_Trend_Sen completed ..."

!================================================================================
Do ivar=1,18
!Do it=1, Ntime
   Do ist=1, Nst
!      print "(40F5.1)", (Events(ivar, ilon, jlat, 15), ilon=1,Nlon)
   enddo
 !  print*,ivar; pause 
Enddo
!================================================================================
CONTAINS
!--------------------------------------------------------------------------------
Subroutine Cal_Trend_Sen
Integer :: ist, it, ireg, N1, Ntime_1, Numyear
Integer, Parameter :: Nidx = 6 
Real :: SenSl, prob
Real, Allocatable :: Wk(:), Xt(:), Xstn(:)   ! Chi tinh tu nam thu 2 tro di (bo nam dau tien)

Ntime_1 = Ntime-Nmon
Allocate (SenIDX(Nst,Nidx), ProbIDX(Nst,Nidx), Wk(Ntime_1), Xt(Nidx), Xstn(Nst))
SenIDX(:,:) = Miss; ProbIDX(:,:) = Miss; NumYear = 10

   Do ist=1, Nst
         Wk(1:Ntime_1)=  PDSI(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,1) = SenSl; ProbIDX(ist,1) = prob
         Wk(1:Ntime_1)=scPDSI(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,2) = SenSl; ProbIDX(ist,2) = prob
         Wk(1:Ntime_1)= SPI01(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,3) = SenSl; ProbIDX(ist,3) = prob
         Wk(1:Ntime_1)= SPI03(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,4) = SenSl; ProbIDX(ist,4) = prob
         Wk(1:Ntime_1)= SPI06(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,5) = SenSl; ProbIDX(ist,5) = prob
         Wk(1:Ntime_1)= SPI12(ist,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ist,6) = SenSl; ProbIDX(ist,6) = prob
         print*, "Station: ", ist
   Enddo
FName = "Sen_Prob_idx_Stn.nc"
Do it=1, Nidx; Xt(it) = it; Enddo;  Do ist=1, Nst; Xstn(ist)=ist; Enddo
Call wr_ncdf_2d (trim(Fname), XStn, Xt, Nst, Nidx, SenIDX, "SenSL", "SenSL", "-", "Station", "Dr_idx", Miss)
Call wr_ncdf_2d_apd (trim(Fname), ProbIDX, "Prob", Nst, Nidx, "Prob", "-", "Station", "Dr_idx", Miss)
End Subroutine Cal_Trend_Sen
!--------------------------------------------------------------------------------
Real Function Ave_Spell (X, N)
Integer :: N, N1, i
Real :: X(N), tmp
N1=0; tmp=0
Do i=1, N
   if (X(i) > 0) then
      N1=N1+X(i); tmp=tmp+X(i)*I
   endif
Enddo
if (N1>0 .AND. N1<N) then; tmp=tmp/real(N1); else; tmp=Miss; endif
Ave_Spell=tmp
End Function Ave_Spell
!--------------------------------------------------------------------------------
Subroutine Drought_Events   ! Xac dinh cac thang Drought (0/1) theo cac nguong, cac IDX
! Event_Dr
Integer :: it, ist, ing, iv, imon, imon1, idx, icf, icm
Integer, Parameter :: N_Nguong = 3, Nvars = 18 ! 18 = 6 chi so x 3 nguong
Real, Parameter :: Nguong(N_Nguong) = (/-1.0, -1.5, -2.0/)
Character (Len=255) :: FName
Character (Len=4), Parameter :: VarName(N_Nguong) = (/"Mod_","Sev_","Ext_"/)
Character (Len=30) :: Varn(Nvars)
Real :: Xvar(Nvars), Xmon(Nmon)
Logical :: Append
Real, Allocatable, Dimension (:,:,:) :: Freq, Mag, Events_1   ! Drought Frequency & Magnitude
Real, Allocatable :: XStn(:)
Real, Allocatable, Dimension (:,:) :: Freq_all, Mag_all     ! Drought Frequency & Magnitude for entire period
Real, Allocatable :: X1(:), M1(:), Ev(:,:)    ! Working array
Real, Allocatable :: XX(:,:)    ! Working array
Real :: xmag, Imax1    ! Do dai lon nhat mot dot han
Integer :: Imax
!Real, Allocatable :: X1(:), M1(:)    ! Working array
Real, Allocatable :: X_Spell(:,:,:), X_Time_Arr(:,:,:)
!Real :: Imax1    ! Do dai lon nhat mot dot han
!Integer :: Imax
!        ----------------------------------------------------
!Varn(1)="pdsi_mod"; Varn(2)="scpdsi_mod"; Varn(3)="spi01_mod"; Varn(4)="spi03_mod"; Varn(5)="spi06_mod"; Varn(6)="spi12_mod"
!Varn(7)="pdsi_sev"; Varn(8)="scpdsi_sev"; Varn(9)="spi01_sev"; Varn(10)="spi03_sev"; Varn(11)="spi06_sev"; Varn(12)="spi12_sev"
!Varn(13)="pdsi_ext"; Varn(14)="scpdsi_ext"; Varn(15)="spi01_ext"; Varn(16)="spi03_ext"; Varn(17)="spi06_ext"; Varn(18)="spi12_ext"
Varn(1)="pdsi_mod"; Varn(2)="scpdsi_mod"; Varn(3)="spei06_mod"; Varn(4)="spei12_mod"; Varn(5)="spi06_mod"; Varn(6)="spi12_mod"
Varn(7)="pdsi_sev"; Varn(8)="scpdsi_sev"; Varn(9)="spei06_sev"; Varn(10)="spei12_sev"; Varn(11)="spi06_sev"; Varn(12)="spi12_sev"
Varn(13)="pdsi_ext"; Varn(14)="scpdsi_ext"; Varn(15)="spei06_ext"; Varn(16)="spei12_ext"; Varn(17)="spi06_ext"; Varn(18)="spi12_ext"

Allocate (Event_Dr(Nst, Ntime), Events(Nvars, Nst, Ntime), XStn(Nst), Xtime(Ntime), Value_IDX(Nvars, Nst, Ntime))  ! 18 = 6 chi so x 3 nguong
Allocate (Ev(Nst, Ntime))   ! Luu tam thoi gia tri Dr_IDX
Value_IDX = Miss   ! Value_IDX luu gia tri cac chi so theo nguong

FName = "Drought_Events_Stn.nc"
Do ist=1, Nst; XStn(ist)=ist; Enddo ; Do it=1, Ntime; Xtime(it)=it; Enddo   ! Mang T gia tao
iv = 0
Do ing=1, N_Nguong
   Event_Dr(:,:) = Miss
   Where (PDSI/=Miss .AND. PDSI<=Nguong(ing)*2); Event_Dr = 1; Else Where (PDSI/=Miss .AND. PDSI>Nguong(ing)*2); Event_Dr = 0; End Where
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"pdsi"
   Ev=Miss; Where (Event_Dr==1); Ev=PDSI; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)

   Event_Dr(:,:) = Miss
   Where (scPDSI/=Miss .AND. scPDSI<=Nguong(ing)*2); Event_Dr = 1; Else Where (scPDSI/=Miss .AND. scPDSI>Nguong(ing)*2); Event_Dr = 0; End Where 
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"scpdsi"
   Ev=Miss; Where (Event_Dr==1); Ev=scPDSI; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)

   Event_Dr(:,:) = Miss
   Where (SPI01/=Miss .AND. SPI01<=Nguong(ing)); Event_Dr = 1; Else Where (SPI01/=Miss .AND. SPI01>Nguong(ing)); Event_Dr = 0; End Where
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"SPI01"
   Ev=Miss; Where (Event_Dr==1); Ev=SPI01; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)

   Event_Dr(:,:) = Miss
   Where (SPI03/=Miss .AND. SPI03<=Nguong(ing)); Event_Dr = 1; Else Where (SPI03/=Miss .AND. SPI03>Nguong(ing)); Event_Dr = 0; End Where
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"SPI03"
   Ev=Miss; Where (Event_Dr==1); Ev=SPI03; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)

   Event_Dr(:,:) = Miss
   Where (SPI06/=Miss .AND. SPI06<=Nguong(ing)); Event_Dr = 1; Else Where (SPI06/=Miss .AND. SPI06>Nguong(ing)); Event_Dr = 0; End Where
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"SPI06"
   Ev=Miss; Where (Event_Dr==1); Ev=SPI06; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)

   Event_Dr(:,:) = Miss
   Where (SPI12/=Miss .AND. SPI12<=Nguong(ing)); Event_Dr = 1; Else Where (SPI12/=Miss .AND. SPI12>Nguong(ing)); Event_Dr = 0; End Where
   iv = iv+1; Events(iv,:,:) = Event_Dr(:,:); Varn(iv) = trim(VarName(ing))//"SPI12"
   Ev=Miss; Where (Event_Dr==1); Ev=SPI12; Endwhere; Value_IDX(iv,:,:) = Ev(:,:)
Enddo
Append = .False.
Do iv=1, Nvars
   Event_Dr(:,:) = Events(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(Fname), Event_Dr, trim(Varn(iv)), Nst, Ntime, "Dr_IDX", "-", "Station", "time", Miss)
!      Call wr_ncdf_3d_apd (trim(FName), Event_Dr, trim(Varn(iv)), Nlon, Nlat, Ntime, Miss)
   else
      Call wr_ncdf_2d (trim(Fname), XStn, Xtime, Nst, Ntime, Event_Dr, trim(Varn(iv)), "Dr_IDX", "-", "Station", "time", Miss)
!      Call wr_ncdf_3d (trim(FName), Xlon, Xlat, Xtime, Nlon, Nlat, Ntime, Event_Dr, trim(Varn(iv)), "months since 1966-01-15 00:00:00", Miss)
   endif
   Append = .True.
Enddo

! Do dai trung binh, cuc dai cac dot han; Khoang thoi gian trung binh, cuc dai giua hai dot han:
!Real, Allocatable, Dimension (:,:,:) :: Spell_Ave, Spell_Max, Time_Arr_Ave, Time_Arr_Max
! NSpell_and_Leng (X1, N, M1, Imax1)

Allocate (X1(Ntime), M1(Ntime), Spell_Ave(Nst,Nvars), Spell_Max(Nst,Nvars))
Allocate (X_Spell(Nvars,Nst,200), X_Time_Arr(Nvars,Nst,200)); X_Spell=Miss; X_Time_Arr=Miss

Do iv=1, Nvars
   Do ist=1, Nst
      X1(1:Ntime)=Events(iv,ist,1:Ntime)
      Call NSpell_and_Leng (X1, Ntime, M1, Imax1); Imax = Int(Imax1)
      X_Spell(iv,ist,1:Imax)=M1(1:Imax)
      Spell_Max(ist,iv) = Imax1; Spell_Ave(ist,iv) = Ave_Spell (M1, Imax)
   Enddo
Enddo

! Khoang thoi gian trung binh, cuc dai giua hai dot han (Time arrival)
Allocate (Events_1(Nvars, Nst, Ntime))
Events_1(:,:,:)=Events(:,:,:)  ! Sao luu lai de hoan chuyen giua 0 va 1 cua Events_1
Where (Events_1==0); Events_1=2; End Where; Where (Events_1==1); Events_1=0; End Where
Where (Events_1==2); Events_1=1; End Where
!------------- Print for testing -------------
!open (3, file="Events.txt", status="Unknown")
!Do ist=1, Nst
!   write (3, "(100I3)") (int(Spell_Max(ist,iv)), int((1,ist,it)), ist=1,Nst)
!Enddo
!Close(3)
!---------------------------------------------
Allocate (Time_Arr_Ave(Nst,Nvars), Time_Arr_Max(Nst,Nvars))
Do iv=1, Nvars
!   print*,"iv=",iv
   Do ist=1, Nst
      X1(1:Ntime)=Events_1(iv,ist,1:Ntime)
!      print "(200I2)", int(X1(:))
      Call NSpell_and_Leng (X1, Ntime, M1, Imax1); Imax = Int(Imax1)
      X_Time_Arr(iv,ist,1:Imax)=M1(1:Imax)
!      print "(70I3)", Imax, (int(M1(it)), it=1,Imax)
      Time_Arr_Max(ist,iv) = Imax1; Time_Arr_Ave(ist,iv) = Ave_Spell (M1, Imax)
!    pause
   Enddo
Enddo

!------------- Print for testing -------------
open (3, file="Events.txt", status="Unknown")
Do ist=1, Nst
!   write (3, "(100I3)") (int(Spell_Max(ist,iv)), int(Time_Arr_Max(ist,iv)), iv=1,Nvars)
   write (3, "(6(2I6,2x))") int(Time_Arr_ave(ist,1)), int(Time_Arr_Max(ist,1)), int(Time_Arr_ave(ist,7)), int(Time_Arr_Max(ist,7)),int(Time_Arr_ave(ist,13)), int(Time_Arr_Max(ist,13)) &
                           ,int(Spell_ave(ist,1)), int(Spell_Max(ist,1)), int(Spell_ave(ist,7)), int(Spell_Max(ist,7)),int(Spell_ave(ist,13)), int(Spell_Max(ist,13))
Enddo
!Do it=1, Ntime
!   write (3, "(200I2)") (int(Events(1,ist,it)),int(Events_1(1,ist,it)),int(Events(7,ist,it)),int(Events_1(7,ist,it)), int(Events(13,ist,it)),int(Events_1(13,ist,it)), ist=1,20)
!Enddo
!write(3,*)"===================="
!Do it=1,200
!   write (3,"(120I3)") (int(X_Spell(1,ist,it)),int(X_Time_Arr(1,ist,it)),int(X_Spell(7,ist,it)),int(X_Time_Arr(7,ist,it)),int(X_Spell(13,ist,it)),int(X_Time_Arr(13,ist,it)), ist=1,20)
!Enddo
Close(3)
!---------------------------------------------


Do iv=1, Nvars; Xvar(iv)=iv; Enddo   ! Mang T gia tao
FName = "Spell_Time_Arr_Stn.nc"
Call wr_ncdf_2d     (trim(FName), XStn, Xvar, Nst, Nvars, Spell_Ave, "Spell_Ave", "Spell_Ave","-","Station","Dr_idx", Miss)
Call wr_ncdf_2d_apd (trim(FName), Spell_Max, "Spell_Max", Nst, Nvars, "Spell_Max","-", "Station","Dr_idx", Miss)
Call wr_ncdf_2d_apd (trim(FName), Time_Arr_Ave, "Time_Arr_Ave", Nst, Nvars, "Time_Arr_Ave","-", "Station","Dr_idx", Miss)
Call wr_ncdf_2d_apd (trim(FName), Time_Arr_Max, "Time_Arr_Max", Nst, Nvars, "Time_Arr_Max","-", "Station","Dr_idx", Miss)
!wr_ncdf_2d (filenm,lon,lat,nlon,nlat,var,var_nm, var_long_name, var_unit, lon_nm, lat_nm, flag)
!wr_ncdf_2d_apd (filenm,var,var_nm,nlon,nlat, var_long_name, var_unit, lon_nm, lat_nm, flag)

Deallocate (Events_1, X1, M1)

! Tinh so thang han trong toan chuoi thoi gian
Deallocate(Event_Dr); Allocate(Event_Dr(Nst, Nvars))
Event_Dr(:,:) = Miss; 
Do iv=1, Nvars
   Do ist=1, Nst
   !   if (Reg(ilon,jlat)>0) then
          Event_Dr(ist,iv) = 0
          Do it=1, Ntime
             if (Events(iv,ist,it)>0) Event_Dr(ist,iv) = Event_Dr(ist,iv)+1
          Enddo
   !   endif 
   Enddo   
Enddo
Do iv=1, Nvars; Xvar(iv)=iv; Enddo   ! Mang T gia tao
FName = "Num_Drou_Mon_Stn.nc"
Call wr_ncdf_2d (trim(Fname), XStn, Xvar, Nst, Nvars, Event_Dr, "N_Drou_mon", "Dr_IDX", "-", "Station", "Index", Miss)
!Call wr_ncdf_3d (trim(FName), Xlon, Xlat, Xvar, Nlon, Nlat, Nvars, Event_Dr, "N_Drou_mon", "months since 1966-01-15 00:00:00", Miss)

! Tinh Fre & Mag cho toan chuoi theo Severity (Freq_all, Mag_all)
Allocate (Freq_all(Nst,Nvars),Mag_all(Nst,Nvars)); Freq_all(:,:)=Miss; Mag_all(:,:)=Miss
Do iv=1, Nvars
   Do ist=1, Nst
      icf=0; xmag=0
      Do it=Nmon+1, Ntime       ! Chi tinh tu nam thu 2 tro di
         if (Events(iv,ist,it)>0) then
            !Freq_all(iv,ilon,jlat)=Freq_all(iv,ilon,jlat)+1.
            icf=icf+1
            xmag=xmag + Value_IDX(iv,ist,it)
         endif
      Enddo
      if (icf>0) then
          Mag_all(ist,iv) = xmag/real(icf); Freq_all(ist,iv)= real(icf)/Real(Ntime-Nmon)
      else
          Mag_all(ist,iv) = 0; Freq_all(ist,iv)= 0; 
      endif
   Enddo
Enddo
FName="Fre_Mag_Stn_all_period.nc";
Do iv=1, Nvars; Xvar(iv)=iv; Enddo   ! Mang T gia tao
Call wr_ncdf_2d     (trim(FName), XStn, Xvar, Nst, Nvars, Freq_all, "Freq", "Freq","-","Station","Dr_idx", Miss)
Call wr_ncdf_2d_apd (trim(FName), Mag_all, "Mag", Nst, Nvars, "Mag","-", "Station","Dr_idx", Miss)
!wr_ncdf_2d (filenm,lon,lat,nlon,nlat,var,var_nm, var_long_name, var_unit, lon_nm, lat_nm, flag)

! Tinh tan suat han trong tung thang theo Severity
Allocate (Freq(Nvars,Nst,Nmon),Mag(Nvars,Nst,Nmon)); Freq(:,:,:)=Miss; Mag(:,:,:)=Miss
Do imon=1, Nmon; Xmon(imon)=imon; Enddo   ! Mang T gia tao
Do iv=1, Nvars
   Do ist=1, Nst
      Freq(iv,ist,:)=0; Mag(iv,ist,:)=0
      Do it=Nmon+1, Ntime       ! Chi tinh tu nam thu 2 tro di
         imon=mod(it,Nmon); if (imon==0) imon=Nmon
         if (Events(iv,ist,it)>0) then
            Freq(iv,ist,imon)=Freq(iv,ist,imon)+1; Mag(iv,ist,imon)=Mag(iv,ist,imon)+Value_IDX(iv,ist,it)
!               idx=mod(iv,6)
!               if      (idx==1) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+  PDSI(ist,it);
!               else if (idx==2) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+scPDSI(ist,it);
!               else if (idx==3) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+ SPI01(ist,it);
!               else if (idx==4) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+ SPI03(ist,it);
!               else if (idx==5) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+ SPI06(ist,it);
!               else if (idx==0) then; Mag(iv,ist,imon)=Mag(iv,ist,imon)+ SPI12(ist,it); endif
         endif
      Enddo
      if (Freq(iv,ist,imon)>0) then
         Mag(iv,ist,imon)=Mag(iv,ist,imon) /Freq(iv,ist,imon)         ! Trung binh SPI/PDSI nhung thang han
         Freq(iv,ist,imon)=Freq(iv,ist,imon)/Real((Ntime-Nmon)/Nmon)  !*100.    ! Tan suan han (% so thang han/Tong so nam)
      else
         Mag(iv,ist,imon)=0; Freq(iv,ist,imon)=0
      endif
   Enddo    
Enddo
FName="Fre_Mag_Stn.nc"
Deallocate(Event_Dr); Allocate(Event_Dr(Nst, Nmon))
Event_Dr(:,:) = Miss;
Append = .False.
Do iv=1, Nvars   ! write Frequency
   Event_Dr(:,:) = Freq(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(Fname), Event_Dr, "Fre_"//trim(Varn(iv)), Nst, Nmon, "Dr_IDX", "-", "Station", "time", Miss)
   else
      Call wr_ncdf_2d (trim(Fname), XStn, Xmon, Nst, Nmon, Event_Dr, "Fre_"//trim(Varn(iv)), "Dr_IDX", "-", "Station", "time", Miss)
   endif
   Append = .True.
Enddo
Do iv=1, Nvars   ! write Magnitude
   Event_Dr(:,:) = Mag(iv,:,:)
   Call wr_ncdf_2d_apd (trim(Fname), Event_Dr, "Mag_"//trim(Varn(iv)), Nst, Nmon, "Dr_IDX", "-", "Station", "time", Miss)
Enddo

End Subroutine Drought_Events
!--------------------------------------------------------------------------------
Subroutine Ave_subreg_PDSI_SPI   ! Trung binh vung KH cac chi so va Phan tram dien tich Drought
Integer :: it, ireg, i, iTh, idx1, idx2, Nidx
Real :: tmp, Nguong
Real :: Xpar(Npar), Ypar(1), Ytime(Ntime) 
Real :: XArea (Nidx_Area)
Real :: XXArea(Nidx_Area/3)   ! De ghi tach thanh 3 file Moderate, Severe, Extreme
Character (Len=2),Parameter :: VG(0:Nreg) = (/"VN","B1","B2","B3","B4","N1","N2","N3"/)
Integer, Parameter :: NThre = 3
Real, Parameter :: Thre(NThre) = (/-1.0, -1.5, -2.0/)
Character (Len=5) :: Severity(3) = (/"Moder","Sever","Extre"/)
Real, Allocatable :: Per_tmp(:,:)

Allocate (Ave_IDX(Npar, Ntime)); Ave_IDX(:,:) = Miss

! Trung binh vung theo cac vung KH cua cac chi so PDSI, scPDSI, SPI01-12 (Npar=42 chi so) 
! (1-7: PDSI, 2-14:scPDSI,...
Do it=1, Ntime
   Do ireg=1, Nreg
      Ave_IDX(ireg,it)=Area_Mean(PDSI,ireg,it);         Ave_IDX(ireg+Nreg,it)  =Area_Mean(scPDSI,ireg,it)
      Ave_IDX(ireg+Nreg*2,it)=Area_Mean(SPI01,ireg,it); Ave_IDX(ireg+Nreg*3,it)=Area_Mean(SPI03,ireg,it)
      Ave_IDX(ireg+Nreg*4,it)=Area_Mean(SPI06,ireg,it); Ave_IDX(ireg+Nreg*5,it)=Area_Mean(SPI12,ireg,it)
   Enddo
Enddo
Do i=1, Npar; Xpar(i)=i; Enddo; Ypar(1)=1
Do i=1, Ntime; Ytime(i)=i; Enddo

FName="Ave_Drought_IDX_Stn.nc"
Call wr_ncdf_2d (trim(Fname), Xpar, Ytime, Npar, Ntime, Ave_IDX, "Ave_Idx", "Dr_IDX", "-", "Index", "time", Miss)
!Call wr_ncdf_3d  (trim(FName), Ypar, Xpar, Ytime, 1, Npar, Ntime, Ave_IDX, "Ave_Idx", "months since 1966-01-15 00:00:00", Miss)

! Phan tram dien tich drought theo Severity: Moderate, Severe, Extreme
! NThre -> Nreg(=8) -> Nidx/3(=6)
Allocate (Percent_Area(Nidx_Area,Ntime)); 

Do it=1, Ntime 
   i=0
   Do iTh=1, NThre  ! ireg=0, Nreg 
      Do ireg=0, Nreg  !iTh=1, NThre
         i=i+1
         Percent_Area(i,  it)= Percent_Dr(PDSI,   it,Thre(iTh)*2, VG(ireg))
         Percent_Area(i+1,it)= Percent_Dr(scPDSI, it,Thre(iTh)*2, VG(ireg))
         Percent_Area(i+2,it)= Percent_Dr(SPI01,  it,Thre(iTh),   VG(ireg))
         Percent_Area(i+3,it)= Percent_Dr(SPI03,  it,Thre(iTh),   VG(ireg))
         Percent_Area(i+4,it)= Percent_Dr(SPI06,  it,Thre(iTh),   VG(ireg))
         Percent_Area(i+5,it)= Percent_Dr(SPI12,  it,Thre(iTh),   VG(ireg))
         i=i+5
      Enddo
   Enddo
  ! print"(18F7.2)",Percent_Area(1,:,it)
Enddo
Nidx = Nidx_Area/3
Do i=1, Nidx_Area; XArea(i)=i; Enddo; Do i=1, Nidx; XXArea(i)=i; Enddo
Allocate (Per_tmp(Nidx,Ntime))
Do i=1,3
   idx1=(i-1)*Nidx+1; idx2=i*Nidx
   Per_tmp(1:Nidx,1:Ntime) = Percent_Area(idx1:idx2,1:Ntime)
   FName="Percent_Dr_Area_"//trim(Severity(i))//"_Stn.nc"
   Call wr_ncdf_2d (trim(Fname), XXArea, Ytime, Nidx, Ntime, Per_tmp, "Dr_Area", "Dr_IDX", "-", "Index", "time", Miss)
!   Call wr_ncdf_3d  (trim(FName), Ypar, XXArea, Ytime, 1, Nidx, Ntime, Per_tmp, "Dr_Area", "months since 1966-01-15 00:00:00", Miss)
Enddo
!FName="Percent_Dr_Area.nc"
!Call wr_ncdf_3d  (trim(FName), Ypar, XArea, Ytime, 1, Nidx_Area, Ntime, Percent_Area, "Dr_Area", "months since 1966-01-15 00:00:00", Miss)

End Subroutine Ave_subreg_PDSI_SPI
!--------------------------------------------------------------------------------
Real Function Percent_Dr (X,it,Nguong,Vung)  ! Phan tram dien tich co drought
Integer :: it, ist, ireg, icount, iDr
Character (Len=*) :: Vung
Real, Dimension(Nst,Ntime) :: X
Real :: tmp, Nguong
Logical :: LVung

tmp = 0; icount=0; iDr=0
Do ist=1, Nst
   
      ireg = INT(Reg(ist))
      if (Vung=="VN") then
         LVung = ireg>0
      else if (Vung=="B1") then
         LVung = ireg==1
      else if (Vung=="B2") then
         LVung = ireg==2
      else if (Vung=="B3") then
         LVung = ireg==3
      else if (Vung=="B4") then
         LVung = ireg==4
      else if (Vung=="N1") then
         LVung = ireg==5
      else if (Vung=="N2") then
         LVung = ireg==6
      else if (Vung=="N3") then
         LVung = ireg==7
      endif
      if (LVung) then !Reg(ilon,jlat)>0) then
         icount=icount+1              ! So o luoi VN
         if (X(ist,it)/=Miss .AND. X(ist,it)<Nguong) then  ! Drought event
            iDr=iDr+1
         endif
   !      print*,ireg,icount,iDr
      endif
   
Enddo
if (icount>0) tmp=real(iDr)/real(icount); Percent_Dr = tmp*100

End Function Percent_Dr
!--------------------------------------------------------------------------------
Real Function Area_Mean (X, ireg, it)   ! Tinh trung binh vung chi so X nao do
Integer :: it, ist, ireg, icount
Real, Dimension(Nst,Ntime) :: X
Real :: tmp
tmp = 0; icount=0
Do ist=1, Nst
   if (INT(Reg(ist))==ireg) then
      if (X(ist,it)/=Miss) then
         tmp=tmp+X(ist,it); icount=icount+1
      endif
   endif
Enddo
if (icount>0) tmp=tmp/real(icount); Area_Mean = tmp
End Function Area_Mean
!--------------------------------------------------------------------------------
Subroutine TQ_PDSI_SPI      ! Tuong quan giua cac chi so
Integer :: it, ist
Real, Allocatable, Dimension (:) :: X, Y
Real, Allocatable :: XStn(:)
Real :: Corr, Prob
Character (Len=100) :: FProb

!  Empty_Ar1 (X, N)
Allocate(TQ_PDSI_SPI01(Nst), TQ_PDSI_SPI03(Nst), TQ_PDSI_SPI06(Nst), TQ_PDSI_SPI12(Nst), XStn(Nst))
Allocate(Pr_PDSI_SPI01(Nst), Pr_PDSI_SPI03(Nst), Pr_PDSI_SPI06(Nst), Pr_PDSI_SPI12(Nst))
Allocate (X(Ntime), Y(Ntime))
! SPI - PDSI
TQ_PDSI_SPI01(:)=Miss; TQ_PDSI_SPI03(:)=Miss; TQ_PDSI_SPI06(:)=Miss; TQ_PDSI_SPI12(:)=Miss
Pr_PDSI_SPI01(:)=Miss; Pr_PDSI_SPI03(:)=Miss; Pr_PDSI_SPI06(:)=Miss; Pr_PDSI_SPI12(:)=Miss
Do ist=1, Nst
   X(1:Ntime)=PDSI(ist,1:Ntime)
   Y(1:Ntime)=SPI01(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI01(ist)=Corr; Pr_PDSI_SPI01(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI01(ist)=HSTQ(X,Y,Ntime)
   Y(1:Ntime)=SPI03(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI03(ist)=Corr; Pr_PDSI_SPI03(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI03(ist)=HSTQ(X,Y,Ntime)
   Y(1:Ntime)=SPI06(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI06(ist)=Corr; Pr_PDSI_SPI06(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI06(ist)=HSTQ(X,Y,Ntime)
   Y(1:Ntime)=SPI12(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI12(ist)=Corr; Pr_PDSI_SPI12(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI12(ist)=HSTQ(X,Y,Ntime)
Enddo
FName="Tq_pdsi_spi_Stn.nc"; FProb="Prob_Tq_pdsi_spi_Stn.nc"
Do ist=1, Nst; XStn(ist)=ist; Enddo
Call wr_ncdf_1d     (trim(FName), XStn, Nst, TQ_PDSI_SPI01, "pdsi_spi01", "HSTQ_pdsi_spi01", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI03, "pdsi_spi03", Nst, "HSTQ_pdsi_spi03", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI06, "pdsi_spi06", Nst, "HSTQ_pdsi_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI12, "pdsi_spi12", Nst, "HSTQ_pdsi_spi12", "-", "Station", Miss)

Call wr_ncdf_1d     (trim(FProb), XStn, Nst, Pr_PDSI_SPI01, "pdsi_spi01", "HSTQ_pdsi_spi01", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI03, "pdsi_spi03", Nst, "HSTQ_pdsi_spi03", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI06, "pdsi_spi06", Nst, "HSTQ_pdsi_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI12, "pdsi_spi12", Nst, "HSTQ_pdsi_spi12", "-", "Station", Miss)


!Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, TQ_PDSI_SPI01, "pdsi_spi01", Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI03, "pdsi_spi03", Nlon,Nlat, Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI06, "pdsi_spi06", Nlon,Nlat, Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI12, "pdsi_spi12", Nlon,Nlat, Miss)

! SPI - scPDSI
TQ_PDSI_SPI01(:)=Miss; TQ_PDSI_SPI03(:)=Miss; TQ_PDSI_SPI06(:)=Miss; TQ_PDSI_SPI12(:)=Miss
Pr_PDSI_SPI01(:)=Miss; Pr_PDSI_SPI03(:)=Miss; Pr_PDSI_SPI06(:)=Miss; Pr_PDSI_SPI12(:)=Miss
Do ist=1, Nst
   X(1:Ntime)=scPDSI(ist,1:Ntime)
   Y(1:Ntime)=SPI01(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI01(ist)=Corr; Pr_PDSI_SPI01(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI01(ist)=HSTQ(X,Y,Ntime)
   Y(1:Ntime)=SPI03(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI03(ist)=Corr; Pr_PDSI_SPI03(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI03(ist)=HSTQ(X,Y,Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI06(ist)=Corr; Pr_PDSI_SPI06(ist)=Prob
   Y(1:Ntime)=SPI06(ist,1:Ntime)
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI06(ist)=HSTQ(X,Y,Ntime)
   Y(1:Ntime)=SPI12(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI12(ist)=Corr; Pr_PDSI_SPI12(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI12(ist)=HSTQ(X,Y,Ntime)
Enddo
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI01, "scpdsi_spi01", Nst, "HSTQ_pdsi_spi01", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI03, "scpdsi_spi03", Nst, "HSTQ_pdsi_spi03", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI06, "scpdsi_spi06", Nst, "HSTQ_pdsi_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI12, "scpdsi_spi12", Nst, "HSTQ_pdsi_spi12", "-", "Station", Miss)

Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI01, "scpdsi_spi01", Nst, "HSTQ_pdsi_spi01", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI03, "scpdsi_spi03", Nst, "HSTQ_pdsi_spi03", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI06, "scpdsi_spi06", Nst, "HSTQ_pdsi_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI12, "scpdsi_spi12", Nst, "HSTQ_pdsi_spi12", "-", "Station", Miss)

! SEPI - SPI
TQ_PDSI_SPI01(:)=Miss; TQ_PDSI_SPI03(:)=Miss; TQ_PDSI_SPI06(:)=Miss; TQ_PDSI_SPI12(:)=Miss
Pr_PDSI_SPI01(:)=Miss; Pr_PDSI_SPI03(:)=Miss; Pr_PDSI_SPI06(:)=Miss; Pr_PDSI_SPI12(:)=Miss
Do ist=1, Nst
   X(1:Ntime)=SPI01(ist,1:Ntime)
   Y(1:Ntime)=SPI06(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI01(ist)=Corr; Pr_PDSI_SPI01(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI01(ist)=HSTQ(X,Y,Ntime)  ! SPEI06 - SPI06
   Y(1:Ntime)=SPI12(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI03(ist)=Corr; Pr_PDSI_SPI03(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI03(ist)=HSTQ(X,Y,Ntime)  ! SPEI06 - SPI12
   X(1:Ntime)=SPI03(ist,1:Ntime)
   Y(1:Ntime)=SPI06(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI06(ist)=Corr; Pr_PDSI_SPI06(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI06(ist)=HSTQ(X,Y,Ntime)  ! SPEI12 - SPI06
   Y(1:Ntime)=SPI12(ist,1:Ntime)
   Call Corre (X, Y, Ntime, Corr, Prob, Miss); TQ_PDSI_SPI12(ist)=Corr; Pr_PDSI_SPI12(ist)=Prob
!   if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI12(ist)=HSTQ(X,Y,Ntime)  ! SPEI12 - SPI12
Enddo
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI01, "spei06_spi06", Nst, "HSTQ_spei06_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI03, "spei06_spi12", Nst, "HSTQ_spei06_spi12", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI06, "spei12_spi06", Nst, "HSTQ_spei12_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI12, "spei12_spi12", Nst, "HSTQ_spei12_spi12", "-", "Station", Miss)

Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI01, "spei06_spi06", Nst, "HSTQ_spei06_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI03, "spei06_spi12", Nst, "HSTQ_spei06_spi12", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI06, "spei12_spi06", Nst, "HSTQ_spei12_spi06", "-", "Station", Miss)
Call wr_ncdf_1d_apd (trim(FProb), Pr_PDSI_SPI12, "spei12_spi12", Nst, "HSTQ_spei12_spi12", "-", "Station", Miss)

!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI01, "scpdsi_spi01", Nlon,Nlat, Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI03, "scpdsi_spi03", Nlon,Nlat, Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI06, "scpdsi_spi06", Nlon,Nlat, Miss)
!Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI12, "scpdsi_spi12", Nlon,Nlat, Miss)

! SPI - Precip   <<<<<<<<<<<<< Chua thuc hien
!Allocate (TQ_Precip_SPI01(Nst), TQ_Precip_SPI03(Nst), TQ_Precip_SPI06(Nst), TQ_Precip_SPI12(Nst))
!TQ_Precip_SPI01(:)=Miss; TQ_Precip_SPI03(:)=Miss; TQ_Precip_SPI06(:)=Miss; TQ_Precip_SPI12(:)=Miss
!Do ist=1, Nst 
!      X(1:Ntime)=Precip(ist,1:Ntime)
!      Y(1:Ntime)=SPI01(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI01(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=SPI03(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI03(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=SPI06(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI06(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=SPI12(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI12(ist)=HSTQ(X,Y,Ntime)
!Enddo
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI01, "scpdsi_spi01", Nst, "HSTQ_pdsi_spi01", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI03, "scpdsi_spi03", Nst, "HSTQ_pdsi_spi03", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI06, "scpdsi_spi06", Nst, "HSTQ_pdsi_spi06", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI12, "scpdsi_spi12", Nst, "HSTQ_pdsi_spi12", "-", "Station", Miss)

! PDSI - Precip
!TQ_Precip_SPI01(:)=Miss; TQ_Precip_SPI03(:)=Miss; TQ_Precip_SPI06(:)=Miss; TQ_Precip_SPI12(:)=Miss
!Do ist=1, Nst
!      X(1:Ntime)=Precip(ist,1:Ntime)
!      Y(1:Ntime)=PDSI(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI01(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=scPDSI(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI03(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=ZIND(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI06(ist)=HSTQ(X,Y,Ntime)
!      Y(1:Ntime)=scZIND(ist,1:Ntime)
!      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI12(ist)=HSTQ(X,Y,Ntime)
!Enddo
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI01, "pre_pdsi", Nst, "HSTQ_pre_pdsi", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI03, "pre_scpdsi", Nst, "HSTQ_pre_scpdsi", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI06, "pre_zind", Nst, "HSTQ_pre_zind", "-", "Station", Miss)
!Call wr_ncdf_1d_apd (trim(Fname), TQ_PDSI_SPI12, "pre_sczind", Nst, "HSTQ_pre_sczind", "-", "Station", Miss)
! <<<<<<<<<<<<<<<<<<<<<
End Subroutine TQ_PDSI_SPI
!------------------------------------------------------------------------------
Real Function HSTQ (X, Y, N)
Integer :: N, i, j, N1
Real :: X(N), Y(N), tbx, tby, sx, sy, rxy, tmp
tbx=0.; tby=0.; sx=0.; sy=0.;  rxy=0.; N1=0
Do i=1, N
   if ((x(i)/=Miss).AND.(y(i)/=Miss)) then
      tbx=tbx+x(i); tby=tby+y(i); sx=sx+x(i)*x(i); sy=sy+y(i)*y(i); rxy=rxy+x(i)*y(i)
      N1=N1+1
   endif
Enddo
if (N1>N/2) then
   tbx=tbx/real(n1); tby=tby/real(n1); sx=sx/real(n1); sy=sy/real(n1); rxy=rxy/real(n1)
   sx=sqrt(sx-tbx*tbx); sy=sqrt(sy-tby*tby); rxy=rxy-tbx*tby
   tmp=rxy/(sx*sy)
else
   tmp=Miss
endif
HSTQ = tmp
End Function HSTQ
!------------------------------------------------------------------------------
Subroutine Read_PDSI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: tmplon, tmplat, tmpTime, XStn
Integer :: Nlatmp, Nlotmp, Ntmp
Character (Len=50) :: VarName

!VarName="ori_pdsi"
!  wr_ncdf_2d (trim(Fname), XStn, t, Nst, Ntime, X1, "ori_pdsi", "ori_pdsi", "-", "Station", "Time", amssng)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, PDSI, "ori_pdsi")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, WPLM, "ori_wplm")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, ZIND, "ori_zind")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, scPDSI, "sel_pdsi")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, scWPLM, "sel_wplm")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, scZIND, "sel_zind")
Deallocate (XStn, tmpTime)

! Write out IDX to Text file
!Call Write_IDX_To_TEXT (PDSI, Nst, Ntime, Lon, Lat, Reg, "PDSI", "The PDSI drought index", "PDSI.txt")

End Subroutine Read_PDSI_nc
!--------------------------------------------------------------------------------
Subroutine Read_SPI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: XStn, tmpTime
Integer ::  Nlon1, Nlat1, Nt

!Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI01, "SPI_01")
!Deallocate (XStn, tmpTime)
!Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI03, "SPI_03")
!Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI06, "SPI_06")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI12, "SPI_12")
Deallocate (XStn, tmpTime)
!print*,SPI12; pause
! rd_ncdf_1d (filenm, lon, nlon, var, var_nm)
Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlon, "lon")
Deallocate (XStn)
Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlat, "lat")
Deallocate (XStn)
Call rd_ncdf_1d (trim(FName), XStn, Nst, Reg,  "reg")
End Subroutine Read_SPI_nc
!--------------------------------------------------------------------------------
Subroutine Read_SPEI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: XStn, tmpTime
Integer ::  Nlon1, Nlat1, Nt

Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI01, "SPI_06")
Deallocate (XStn, tmpTime)
Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI03, "SPI_12")
Deallocate (XStn, tmpTime)
!Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI06, "SPI_06")
!Deallocate (XStn, tmpTime)
!Call rd_ncdf_2d (trim(FName), XStn, tmpTime, Nst, Ntime, SPI12, "SPI_12")
!Deallocate (XStn, tmpTime)
!print*,SPI12; pause
! rd_ncdf_1d (filenm, lon, nlon, var, var_nm)
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlon, "lon")
!Deallocate (XStn)
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Xlat, "lat")
!Deallocate (XStn)
!Call rd_ncdf_1d (trim(FName), XStn, Nst, Reg,  "reg")

! wr_ncdf_1d_apd (trim(Fname), lon, "lon", Nst, "longitude", "deg. E", "Station", amssng)

End Subroutine Read_SPEI_nc
!------------------------------------------------------------------------------
Logical Function Empty_Ar1 (X, N)
! Tra ve .True. neu X(N) khuyet toan bo SL
Integer :: N, i, k
Real :: X(N)
k = 0 ! Dem so phan tu khuyet cua X
Do i=1, N
   if (X(i) == Miss) k = k + 1
Enddo
if (k > N/10) then ; Empty_Ar1 = .True. ; else ; Empty_Ar1 = .False. ; endif
End Function Empty_Ar1

!------------------------------------------------------------------------------
Subroutine Corre (X, Y, N, Corr, Prob, Miss)
Integer :: N, i, j, N1
Real :: Corr, Prob, T_Dist
Real :: X(N), Y(N), tbx, tby, sx, sy, rxy, Miss, t
tbx=0.; tby=0.; sx=0.; sy=0.;  rxy=0.; n1=0
Do i=1, N
   if (x(i)/=Miss.And.y(i)/=Miss) then
      n1=n1+1
      tbx=tbx+x(i); tby=tby+y(i); sx=sx+x(i)*x(i); sy=sy+y(i)*y(i); rxy=rxy+x(i)*y(i)
   endif
Enddo
tbx=tbx/real(n1); tby=tby/real(n1); sx=sx/real(n1); sy=sy/real(n1); rxy=rxy/real(n1)
sx=sqrt(sx-tbx*tbx); sy=sqrt(sy-tby*tby); rxy=rxy-tbx*tby
Corr = rxy/(sx*sy)
if (N1<=2) then
   Corr = Miss; Prob = Miss
else
   t = ABS(Corr/ ( Sqrt(1-Corr**2) / Sqrt(Real(N1-2)) ) )
   Prob = T_Dist (t, Real(N1-2))
!   if (Prob<0.OR.Prob>0.5) print*,"T_Dist:", t, N1-2, Prob
endif
End Subroutine Corre
!--------------------------------------------------------------------------------
!---------------------------------------------------------------------!
Subroutine Write_IDX_To_TEXT (X, Nst, Ntime, Lon, Lat, Reg, Varname, Title, Fname)
Integer, Parameter :: Yr1=1980, Yr2=2014
Integer :: Nst, Ntime, ist, iyr, imon, it, ireg, Nyear
Character (Len=*) :: Varname, Title, Fname
Character (Len=50) :: Fmt
Real :: X(Nst, Ntime), Lon(Nst), Lat(Nst), Reg(Nst)

Open (3, file=trim(Fname), status="Unknown")
write (Fmt, "(A,I3,A)") "(I4.4,I2.2,",Nst,"F8.3)"  ! yyyymm, Nst values
write (3,"(A)") Title
write (3, "(A6,200(I6,2x))") "Reg=", (Int(Reg(ist)), ist=1, Nst)
write (3, "(A6,200F8.3)") "Lon=", (Lon(ist), ist=1, Nst)
write (3, "(A6,200F8.3)") "Lat=", (Lat(ist), ist=1, Nst)
Nyear=Ntime/Nmon;
if (Nyear /= (Yr2-Yr1+1)) then; Print*, "Loi !!! Kiem tra lai chuoi so lieu ..."; return; endif  
it=0
Do iyr=Yr1, Yr2
   Do imon=1, Nmon
      it=it+1
      write (3,Fmt) iyr, imon, (X(ist, it), ist=1, Nst)
   Enddo
Enddo
Close(3)
End Subroutine Write_IDX_To_TEXT
!---------------------------------------------------------------------!

Include "../Libs/sub_rd_ncdf_Tan.F90"
!Include "../Libs/sub_rd_ncdf.F90"
Include "../Libs/sub_wr_ncdf_Tan.F90"
Include "../Sub_Trend_MK_Linear.f90"
Include "Num_Dr_Spells.f90"
END 
!----------------------------------------------------------------------------------!
Include "/home/tanpv/SEACLID/Libs/Student_Prob.f90"
!Include "ENSO_FA.f90"
!----------------------------------------------------------------------------------!

!     ori_pdsi  0  t,y,x  What the hell (!@.@) !!!
!     ori_wplm  0  t,y,x  What the hell (!@.@) !!!
!     ori_zind  0  t,y,x  What the hell (!@.@) !!!
!     sel_pdsi  0  t,y,x  What the hell (!@.@) !!!
!     sel_wplm  0  t,y,x  What the hell (!@.@) !!!
!     sel_zind  0  t,y,x  What the hell (!@.@) !!!
!     subreg  0  y,x  What the hell (!@.@) !!!