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
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
! 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/RES/PDSI_Res/PDSI_Stn.nc"
   Call Read_PDSI_nc (Fname)
! Reading SPI
   Fname="/home/tanpv/SEACLID/RES/SPI_Res/SPI_Stn.nc"
   Call Read_SPI_nc (FName)
!================================================================================
   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 :: X1(:), M1(:)    ! Working array
Real, Allocatable :: X_Spell(:,:,:), X_Time_Arr(:,:,:)
Real :: Imax1    ! Do dai lon nhat mot dot han
Integer :: Imax
!        ----------------------------------------------------

Allocate (Event_Dr(Nst, Ntime), Events(Nvars, Nst, Ntime), XStn(Nst), Xtime(Ntime))  ! 18 = 6 chi so x 3 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"
   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"
   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"
   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"
   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"
   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"
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 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
   !   if (Reg(ilon,jlat)>0) then
         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
               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
         Mag(iv,ist,imon)=Mag(iv,ist,imon)  !/Freq(iv,ilon,jlat,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)
   !   endif  
   Enddo    
Enddo
FName="Fre_Mag_Stn.nc"
Deallocate(Event_Dr); Allocate(Event_Dr(Nst, Nmon))
Event_Dr(:,:) = Miss;
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"

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)
!      Call wr_ncdf_3d_apd (trim(FName), Event_Dr, "Fre_"//trim(Varn(iv)), Nlon, Nlat, Nmon, Miss)
   else
      Call wr_ncdf_2d (trim(Fname), XStn, Xmon, Nst, Nmon, Event_Dr, "Fre_"//trim(Varn(iv)), "Dr_IDX", "-", "Station", "time", Miss)
!      Call wr_ncdf_3d (trim(FName), Xlon, Xlat, Xmon, Nlon, Nlat, Nmon, Event_Dr, "Fre_"//trim(Varn(iv)), "months since 1966-01-15 00:00:00", 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)
!   Call wr_ncdf_3d_apd (trim(FName), Event_Dr, "Mag_"//trim(Varn(iv)), Nlon, Nlat, Nmon, 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) 
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
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(:)
!  Empty_Ar1 (X, N)
Allocate(TQ_PDSI_SPI01(Nst), TQ_PDSI_SPI03(Nst), TQ_PDSI_SPI06(Nst), TQ_PDSI_SPI12(Nst), XStn(Nst))
Allocate (X(Ntime), Y(Ntime))
! SPI - PDSI
TQ_PDSI_SPI01(:)=Miss; TQ_PDSI_SPI03(:)=Miss; TQ_PDSI_SPI06(:)=Miss; TQ_PDSI_SPI12(:)=Miss
Do ist=1, Nst
   X(1:Ntime)=PDSI(ist,1:Ntime)
   Y(1:Ntime)=SPI01(ist,1:Ntime)
   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)
   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)
   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)
   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"
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_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
Do ist=1, Nst
   X(1:Ntime)=scPDSI(ist,1:Ntime)
   Y(1:Ntime)=SPI01(ist,1:Ntime)
   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)
   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)
   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)
   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_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)
!print*,PDSI; pause
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")

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

End Subroutine Read_SPI_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
!--------------------------------------------------------------------------------

Include "../Libs/sub_rd_ncdf_Tan.F90"
Include "../Libs/sub_wr_ncdf_Tan.F90"
Include "../Sub_Trend_MK_Linear.f90"
Include "Num_Dr_Spells.f90"
END 

!     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 (!@.@) !!!