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, jlat, ilon, it, Ntime, Nlat, Nlon
!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.nc"
   Call Read_PDSI_nc (Fname)
! Reading SPI
   Fname="/home/tanpv/SEACLID/RES/SPI_Res/SPEI.nc"
   Call Read_SPI_nc (FName)
! Reading Vietnam Mask
   Fname="/home/tanpv/SEACLID/MASK/VN.grid_0.25.nc"
   Call Read_VN_nc (FName)
!================================================================================
   Call TQ_PDSI_SPI            ! Write out "Tq_pdsi_spi.nc"
   Call Ave_subreg_PDSI_SPI    ! Write out "Ave_Drought_IDX.nc", "Percent_Dr_Area.nc" 
   Call Drought_Events         ! Write out "Drought_Events.nc", "Num_Drou_Mon.nc", "Fre_Mag.nc"
   Call Cal_Trend_Sen

!================================================================================
Do ivar=1,18
!Do it=1, Ntime
   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (Events(ivar, ilon, jlat, 15), ilon=1,Nlon)
   enddo
 !  print*,ivar; pause 
Enddo
   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (TQ_PDSI_SPI03(ilon,jlat), ilon=1,Nlon)
   enddo
!   pause
   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (TQ_PDSI_SPI06(ilon,jlat), ilon=1,Nlon)
   enddo
!   pause
   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (TQ_PDSI_SPI12(ilon,jlat), ilon=1,Nlon)
   enddo
!      print*,"PDSI, it=",it
!   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (Precip(ilon,jlat,it), ilon=1,Nlon)
!   Enddo
!      print*,"SPI, it=",it
!   Do jlat=Nlat,1, -1
!      print "(40F5.1)", (Reg(ilon,jlat), ilon=1,Nlon)
!   Enddo
!      print*,"VN - Climatic subregions 1..7"
!   pause
!Enddo
!--------------------------------------------------------------------------------
!================================================================================
CONTAINS
!--------------------------------------------------------------------------------
Subroutine Cal_Trend_Sen
Integer :: jlat, ilon, it, ireg, N1, Ntime_1, Numyear
Integer, Parameter :: Nidx = 6 
Real :: SenSl, prob
Real, Allocatable :: Wk(:), Xt(:)   ! Chi tinh tu nam thu 2 tro di (bo nam dau tien)

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

Do jlat=1, Nlat
   Do ilon=1, Nlon
      if (Reg(ilon,jlat)>0) then
         Wk(1:Ntime_1)=  PDSI(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,1) = SenSl; ProbIDX(ilon,jlat,1) = prob
         Wk(1:Ntime_1)=scPDSI(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,2) = SenSl; ProbIDX(ilon,jlat,2) = prob
         Wk(1:Ntime_1)= SPI01(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,3) = SenSl; ProbIDX(ilon,jlat,3) = prob
         Wk(1:Ntime_1)= SPI03(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,4) = SenSl; ProbIDX(ilon,jlat,4) = prob
         Wk(1:Ntime_1)= SPI06(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,5) = SenSl; ProbIDX(ilon,jlat,5) = prob
         Wk(1:Ntime_1)= SPI12(ilon,jlat,Nmon+1:Ntime); Call MK_Sen(Wk,Ntime_1,NumYear,Miss,SenSl,prob,N1)
         SenIDX(ilon,jlat,6) = SenSl; ProbIDX(ilon,jlat,6) = prob
         print*, "Point at Lon, Lat: ", ilon, jlat
      endif
   Enddo
Enddo
FName = "Sen_Prob_idx.nc"
Do it=1, Nidx; Xt(it) = it; Enddo
Call wr_ncdf_3d (trim(FName), Xlon, Xlat, Xt, Nlon, Nlat, Nidx, SenIDX, "SenSL","months since 1966-01-15 00:00:00","SenSL","-",Miss)
Call wr_ncdf_3d_apd (trim(FName), ProbIDX, "Prob", Nlon, Nlat, Nidx, "Prob", "-", 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) 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, jlat, ilon, 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 :: X1(:), M1(:)    ! Working array
Real, Allocatable :: XX(:,:)    ! Working array
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"

Allocate (Event_Dr(Nlon, Nlat, Ntime), Events(Nvars, Nlon, Nlat, Ntime))  ! 18 = 6 chi so x 3 nguong

FName = "Drought_Events.nc"
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_3d_apd (trim(FName), Event_Dr, trim(Varn(iv)), Nlon, Nlat, Ntime,trim(Varn(iv)),"-", Miss)
   else
      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", trim(Varn(iv)),"-", 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(Nvars, Nlon, Nlat), Spell_Max(Nvars, Nlon, Nlat))
Do iv=1, Nvars
   Do jlat=1, Nlat
      Do ilon=1, Nlon
         if (Reg(ilon,jlat)>0) then
             X1(1:Ntime)=Events(iv,ilon,jlat,1:Ntime)
             Call NSpell_and_Leng (X1, Ntime, M1, Imax1); Imax = Int(Imax1)
             Spell_Max(iv, ilon, jlat) = Imax1; Spell_Ave(iv, ilon, jlat) = Ave_Spell (M1, Imax)
         endif
      Enddo
   Enddo
Enddo

! Khoang thoi gian trung binh, cuc dai giua hai dot han (Time arrival)
Allocate (Events_1(Nvars, Nlon, Nlat, 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

Allocate (Time_Arr_Ave(Nvars, Nlon, Nlat), Time_Arr_Max(Nvars, Nlon, Nlat))
Do iv=1, Nvars
   Do jlat=1, Nlat
      Do ilon=1, Nlon
         if (Reg(ilon,jlat)>0) then
             X1(1:Ntime)=Events_1(iv,ilon,jlat,1:Ntime)
             Call NSpell_and_Leng (X1, Ntime, M1, Imax1); Imax = Int(Imax1)
             Time_Arr_Max(iv, ilon, jlat) = Imax1; Time_Arr_Ave(iv, ilon, jlat) = Ave_Spell (M1, Imax)
         endif
      Enddo
   Enddo
Enddo

FName = "Spell_Ave.nc"; Append = .False.; Allocate(XX(Nlon,Nlat))
Do iv=1, Nvars
   XX(:,:)=Spell_Ave(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(FName), XX, trim(Varn(iv)), Nlon,Nlat,trim(Varn(iv)),"-", Miss)
   else
      Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, XX, trim(Varn(iv)), trim(Varn(iv)),"-", Miss)
      Append = .True.
   endif
Enddo
FName = "Spell_Max.nc"; Append = .False.
Do iv=1, Nvars
   XX(:,:)=Spell_Max(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(FName), XX, trim(Varn(iv)), Nlon,Nlat, trim(Varn(iv)),"-", Miss)
   else
      Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, XX, trim(Varn(iv)), trim(Varn(iv)),"-", Miss)
      Append = .True.
   endif
Enddo
FName = "Time_Arr_Ave.nc"; Append = .False.
Do iv=1, Nvars
   XX(:,:)=Time_Arr_Ave(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(FName), XX, trim(Varn(iv)), Nlon,Nlat, trim(Varn(iv)),"-", Miss)
   else
      Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, XX, trim(Varn(iv)), trim(Varn(iv)),"-", Miss)
      Append = .True.
   endif
Enddo
FName = "Time_Arr_Max.nc"; Append = .False.
Do iv=1, Nvars
   XX(:,:)=Time_Arr_Max(iv,:,:)
   if (Append) then
      Call wr_ncdf_2d_apd (trim(FName), XX, trim(Varn(iv)), Nlon,Nlat, trim(Varn(iv)),"-", Miss)
   else
      Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, XX, trim(Varn(iv)), trim(Varn(iv)),"-", Miss)
      Append = .True.
   endif
Enddo
Deallocate (Events_1, XX, X1, M1)

! Tinh so thang han trong toan chuoi thoi gian
Deallocate(Event_Dr); Allocate(Event_Dr(Nlon, Nlat, Nvars))
Event_Dr(:,:,:) = Miss; 
Do iv=1, Nvars
   Do jlat=1, Nlat
      Do ilon=1, Nlon
         if (Reg(ilon,jlat)>0) then
             Event_Dr(ilon,jlat,iv) = 0
             Do it=1, Ntime
                if (Events(iv,ilon,jlat,it)>0) Event_Dr(ilon,jlat,iv) = Event_Dr(ilon,jlat,iv)+1
             Enddo
         endif 
      Enddo
   Enddo   
Enddo
!
Do iv=1, Nvars; Xvar(iv)=iv; Enddo   ! Mang T gia tao
FName = "Num_Drou_Mon.nc"
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", "N_Drou_mon","-",Miss)

! Tinh tan suat han trong tung thang theo Severity
Allocate (Freq(Nvars,Nlon, Nlat,Nmon),Mag(Nvars,Nlon, Nlat,Nmon)); Freq(:,:,:,:)=Miss; Mag(:,:,:,:)=Miss
Do imon=1, Nmon; Xmon(imon)=imon; Enddo   ! Mang T gia tao
Do iv=1, Nvars
   Do jlat=1, Nlat
      Do ilon=1, Nlon
         if (Reg(ilon,jlat)>0) then
            Freq(iv,ilon,jlat,:)=0; Mag(iv,ilon,jlat,:)=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,ilon,jlat,it)>0) then
                  Freq(iv,ilon,jlat,imon)=Freq(iv,ilon,jlat,imon)+1
                  idx=mod(iv,6)
                  if      (idx==1) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+  PDSI(ilon,jlat,it);
                  else if (idx==2) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+scPDSI(ilon,jlat,it);
                  else if (idx==3) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+ SPI01(ilon,jlat,it);
                  else if (idx==4) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+ SPI03(ilon,jlat,it);
                  else if (idx==5) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+ SPI06(ilon,jlat,it);
                  else if (idx==0) then; Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)+ SPI12(ilon,jlat,it); endif
               endif
            Enddo
            Mag(iv,ilon,jlat,imon)=Mag(iv,ilon,jlat,imon)  !/Freq(iv,ilon,jlat,imon)         ! Trung binh SPI/PDSI nhung thang han
            Freq(iv,ilon,jlat,imon)=Freq(iv,ilon,jlat,imon)!/Real((Ntime-Nmon)/Nmon)*100.    ! Tan suan han (% so thang han/Tong so nam)
         endif
      Enddo
   Enddo    
Enddo
FName="Fre_Mag.nc"
Deallocate(Event_Dr); Allocate(Event_Dr(Nlon, Nlat, 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_3d_apd (trim(FName), Event_Dr, "Fre_"//trim(Varn(iv)), Nlon, Nlat, Nmon, "Fre_"//trim(Varn(iv)), "-",Miss)
   else
      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","Fre_"//trim(Varn(iv)), "-", Miss)
   endif
   Append = .True.
Enddo
Do iv=1, Nvars   ! write Magnitude
   Event_Dr(:,:,:) = Mag(iv,:,:,:)
   Call wr_ncdf_3d_apd (trim(FName), Event_Dr, "Mag_"//trim(Varn(iv)), Nlon, Nlat, Nmon, "Mag_"//trim(Varn(iv)),"-",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(1, Npar, Ntime)); Ave_IDX(1,:,:) = 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(1,ireg,it)=Area_Mean(PDSI,ireg,it);         Ave_IDX(1,ireg+Nreg,it)  =Area_Mean(scPDSI,ireg,it)
      Ave_IDX(1,ireg+Nreg*2,it)=Area_Mean(SPI01,ireg,it); Ave_IDX(1,ireg+Nreg*3,it)=Area_Mean(SPI03,ireg,it)
      Ave_IDX(1,ireg+Nreg*4,it)=Area_Mean(SPI06,ireg,it); Ave_IDX(1,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.nc"
Call wr_ncdf_3d  (trim(FName), Ypar, Xpar, Ytime, 1, Npar, Ntime, Ave_IDX, "Ave_Idx", "months since 1966-01-15 00:00:00", "Ave_Idx","-",Miss)

! Phan tram dien tich drought theo Severity: Moderate, Severe, Extreme
Allocate (Percent_Area(1,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(1,i,  it)= Percent_Dr(PDSI,   it,Thre(iTh)*2, VG(ireg))
         Percent_Area(1,i+1,it)= Percent_Dr(scPDSI, it,Thre(iTh)*2, VG(ireg))
         Percent_Area(1,i+2,it)= Percent_Dr(SPI01,  it,Thre(iTh),   VG(ireg))
         Percent_Area(1,i+3,it)= Percent_Dr(SPI03,  it,Thre(iTh),   VG(ireg))
         Percent_Area(1,i+4,it)= Percent_Dr(SPI06,  it,Thre(iTh),   VG(ireg))
         Percent_Area(1,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(1,Nidx,Ntime))
Do i=1,3
   idx1=(i-1)*Nidx+1; idx2=i*Nidx
   Per_tmp(1,1:Nidx,1:Ntime) = Percent_Area(1,idx1:idx2,1:Ntime)
   FName="Percent_Dr_Area_"//trim(Severity(i))//".nc"
   Call wr_ncdf_3d  (trim(FName), Ypar, XXArea, Ytime, 1, Nidx, Ntime, Per_tmp, "Dr_Area", "months since 1966-01-15 00:00:00", "Dr_Area","-",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, jlat, ilon, ireg, icount, iDr
Character (Len=*) :: Vung
Real, Dimension(Nlon,Nlat,Ntime) :: X
Real :: tmp, Nguong
Logical :: LVung

tmp = 0; icount=0; iDr=0
Do jlat=1, Nlat
   Do ilon=1, Nlon
      ireg = INT(Reg(ilon,jlat))
      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(ilon,jlat,it)/=Miss .AND. X(ilon,jlat,it)<Nguong) then  ! Drought event
            iDr=iDr+1
         endif
   !      print*,ireg,icount,iDr
      endif
   Enddo
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, jlat, ilon, ireg, icount
Real, Dimension(Nlon,Nlat,Ntime) :: X
Real :: tmp
tmp = 0; icount=0
Do jlat=1, Nlat
   Do ilon=1, Nlon
      if (INT(Reg(ilon,jlat))==ireg) then
         if (X(ilon,jlat,it)/=Miss) then
            tmp=tmp+X(ilon,jlat,it); icount=icount+1
         endif
      endif
   Enddo
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, jlat, ilon
Real, Allocatable, Dimension (:) :: X, Y
!  Empty_Ar1 (X, N)
Allocate(TQ_PDSI_SPI01(Nlon,Nlat), TQ_PDSI_SPI03(Nlon,Nlat), TQ_PDSI_SPI06(Nlon,Nlat), TQ_PDSI_SPI12(Nlon,Nlat))
Allocate (X(Ntime), Y(Ntime))
! SPI - PDSI
TQ_PDSI_SPI01(:,:)=Miss; TQ_PDSI_SPI03(:,:)=Miss; TQ_PDSI_SPI06(:,:)=Miss; TQ_PDSI_SPI12(:,:)=Miss
Do jlat=1, Nlat
   Do ilon=1, Nlon
      X(1:Ntime)=PDSI(ilon,jlat,1:Ntime)
      Y(1:Ntime)=SPI01(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI01(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI03(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI03(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI06(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI06(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI12(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI12(ilon,jlat)=HSTQ(X,Y,Ntime)
   Enddo
Enddo
FName="Tq_pdsi_spi.nc"
Call wr_ncdf_2d     (trim(FName), Xlon,Xlat,Nlon,Nlat, TQ_PDSI_SPI01, "pdsi_spi01","pdsi_spi01","-", Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI03, "pdsi_spi03", Nlon,Nlat, "pdsi_spi03", "-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI06, "pdsi_spi06", Nlon,Nlat, "pdsi_spi06","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI12, "pdsi_spi12", Nlon,Nlat, "pdsi_spi12","-",Miss)

! SPI - scPDSI
TQ_PDSI_SPI01(:,:)=Miss; TQ_PDSI_SPI03(:,:)=Miss; TQ_PDSI_SPI06(:,:)=Miss; TQ_PDSI_SPI12(:,:)=Miss
Do jlat=1, Nlat
   Do ilon=1, Nlon
      X(1:Ntime)=scPDSI(ilon,jlat,1:Ntime)
      Y(1:Ntime)=SPI01(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI01(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI03(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI03(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI06(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI06(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI12(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_PDSI_SPI12(ilon,jlat)=HSTQ(X,Y,Ntime)
   Enddo
Enddo
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI01, "scpdsi_spi01", Nlon,Nlat, "scpdsi_spi01","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI03, "scpdsi_spi03", Nlon,Nlat, "scpdsi_spi03","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI06, "scpdsi_spi06", Nlon,Nlat, "scpdsi_spi06","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_PDSI_SPI12, "scpdsi_spi12", Nlon,Nlat, "scpdsi_spi12","-",Miss)

! SPI - Precip
Allocate (TQ_Precip_SPI01(Nlon,Nlat), TQ_Precip_SPI03(Nlon,Nlat), TQ_Precip_SPI06(Nlon,Nlat), TQ_Precip_SPI12(Nlon,Nlat))
TQ_Precip_SPI01(:,:)=Miss; TQ_Precip_SPI03(:,:)=Miss; TQ_Precip_SPI06(:,:)=Miss; TQ_Precip_SPI12(:,:)=Miss
Do jlat=1, Nlat
   Do ilon=1, Nlon
      X(1:Ntime)=Precip(ilon,jlat,1:Ntime)
      Y(1:Ntime)=SPI01(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI01(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI03(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI03(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI06(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI06(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=SPI12(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI12(ilon,jlat)=HSTQ(X,Y,Ntime)
   Enddo
Enddo
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI01, "pre_spi01", Nlon,Nlat, "pre_spi01","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI03, "pre_spi03", Nlon,Nlat, "pre_spi03","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI06, "pre_spi06", Nlon,Nlat, "pre_spi06","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI12, "pre_spi12", Nlon,Nlat, "pre_spi12","-",Miss)

! PDSI - Precip
TQ_Precip_SPI01(:,:)=Miss; TQ_Precip_SPI03(:,:)=Miss; TQ_Precip_SPI06(:,:)=Miss; TQ_Precip_SPI12(:,:)=Miss
Do jlat=1, Nlat
   Do ilon=1, Nlon
      X(1:Ntime)=Precip(ilon,jlat,1:Ntime)
      Y(1:Ntime)=PDSI(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI01(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=scPDSI(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI03(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=ZIND(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI06(ilon,jlat)=HSTQ(X,Y,Ntime)
      Y(1:Ntime)=scZIND(ilon,jlat,1:Ntime)
      if (.not.Empty_Ar1(X,Ntime) .AND. .not.Empty_Ar1(Y,Ntime)) TQ_Precip_SPI12(ilon,jlat)=HSTQ(X,Y,Ntime)
   Enddo
Enddo
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI01, "pre_pdsi", Nlon,Nlat,"pre_pdsi", "-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI03, "pre_scpdsi", Nlon,Nlat, "pre_scpdsi","-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI06, "pre_zind", Nlon,Nlat,"pre_zind", "-",Miss)
Call wr_ncdf_2d_apd (trim(FName), TQ_Precip_SPI12, "pre_sczind", Nlon,Nlat, "pre_sczind","-",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_VN_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: tmplon, tmplat
Integer :: Nlatmp, Nlotmp, Ntmp

Call rd_ncdf_2d (trim(FName), tmplon, tmplat, Nlotmp, Nlatmp, Reg, "all")
Where (Reg==2) Reg=1; Where (Reg==3) Reg=2; Where (Reg==4) Reg=3;
Where (Reg==5) Reg=4; Where (Reg==6) Reg=5; Where (Reg==7) Reg=6;
Where (Reg==8) Reg=7   ! Gan lai cac vung KH: B1=1, B2=2, B3=3, B4=4, N1=5, N2=6, N3=7
Deallocate (tmplon, tmplat)
End Subroutine Read_VN_nc
!--------------------------------------------------------------------------------
Subroutine Read_PDSI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: tmplon, tmplat, tmpTime
Integer :: Nlatmp, Nlotmp, Ntmp
Character (Len=50) :: VarName

VarName="ori_pdsi"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, PDSI, trim(VarName))
Allocate (Xlon(Nlotmp), Xlat(Nlatmp), Xtime(Ntmp))
Xlon(:)=tmplon(:); Xlat(:)=tmplat(:); Xtime(:)=tmpTime(:); Ntime=Ntmp; Nlon=Nlotmp; Nlat=Nlatmp
Deallocate (tmplon, tmplat, tmpTime)
VarName="ori_wplm"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, WPLM, trim(VarName))
Deallocate (tmplon, tmplat, tmpTime)
VarName="ori_zind"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, ZIND, trim(VarName))
Deallocate (tmplon, tmplat, tmpTime)

VarName="sel_pdsi"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scPDSI, trim(VarName))
Deallocate (tmplon, tmplat, tmpTime)
VarName="sel_wplm"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scWPLM, trim(VarName))
Deallocate (tmplon, tmplat, tmpTime)
VarName="sel_zind"
Call rd_ncdf_3d (trim(FName), tmplon, tmplat, tmpTime, Nlotmp, Nlatmp, Ntmp, scZIND, trim(VarName))
Deallocate (tmplon, tmplat, tmpTime)
End Subroutine Read_PDSI_nc
!--------------------------------------------------------------------------------
Subroutine Read_SPI_nc (FName)
Character (Len=*) :: FName
Real, Allocatable, Dimension (:) :: Lon, Lat, Time
Integer ::  Nlon1, Nlat1, Nt

Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI01, "SPI_01")
Xlon(:)=lon(:); Xlat(:)=lat(:); Xtime(:)=Time(:); Ntime=Nt; Nlon=Nlon1; Nlat=Nlat1
Deallocate (Lon, Lat, Time)
Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI03, "SPI_03")
Deallocate (Lon, Lat, Time)
Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI06, "SPI_06")
Deallocate (Lon, Lat, Time)
Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, SPI12, "SPI_12")
Deallocate (Lon, Lat, Time)
Call rd_ncdf_3d (trim(FName), Lon, Lat, Time, Nlon1, Nlat1, Nt, Precip, "Precip")
Deallocate (Lon, Lat, Time)
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.F90"
Include "../Libs/sub_wr_ncdf.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 (!@.@) !!!