! Tap hop cac chuong trinh con tinh cac dac trung Han
!------------------------------------------------------------------------------
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) == amssng) k = k + 1
Enddo
if (k > N/2) then ; Empty_Ar1 = .True. ; else ; Empty_Ar1 = .False. ; endif
End Function Empty_Ar1
!--------------------------------------------------------------
Subroutine Cal_Trend
! MK_Sen(Y,N,Miss,SenSlope,prob,N1)
! Calculate Sen's Slope of time series of:
! 1) Num of Drought months in each year: NDYrs(ist,1:Nyear,ilen)
! 2) Sum of SPI<Nguong in each year: SumSPIYrs(ist,1:Nyear,ilen)
!Real, Dimension(Nst_Max,Nyear,NLen) :: NDYrs ! so thang han cua tung nam ung voi cac qui mo thoi gian
!Real, Dimension(Nst_Max,Nyear,NLen) :: SumSPIYrs ! tong SPI<Nguong cua tung nam ung voi cac qui mo thoi gian
!Real, Dimension(Nst_Max,Nyear,NLen) :: SPI_3m ! tong SPI cua 3 thang nao do trong nam (vi du 2-3-4)
!Real, Dimension(Nst_Max,NLen) :: SenSl_NDYrs,Prob_NDYrs,SenSl_SumSPIYrs,Prob_SumSPIYrs
!Real, Dimension(Nst_Max,NLen) :: SenSl_SPI_3m,Prob_SPI_3m
Integer :: ist, ilen, N1
Real :: wk(Ntime), SenSLp, Prob
NumYear = 20  ! It nhat co 20 nam SL
Do ist=1, Nst
   print*,"Calculating Sen's Slope for point ",Station(ist)
   Do ilen=1, NLen
      wk(1:Nyear-1) = NDYrs(ist,2:Nyear,ilen)
      Call MK_Sen(wk, Nyear-1, NumYear, amssng, SenSlp,Prob,N1)
      SenSl_NDYrs(ist,ilen) = SenSlp; Prob_NDYrs(ist,ilen) = Prob

      wk(1:Nyear-1) = SumSPIYrs(ist,2:Nyear,ilen)
      Call MK_Sen(wk, Nyear-1, NumYear, amssng, SenSlp,Prob,N1)
      SenSl_SumSPIYrs(ist,ilen) = SenSlp; Prob_SumSPIYrs(ist,ilen) = Prob

      wk(1:Nyear-1) = SPI_3m(ist,2:Nyear,ilen)
      Call MK_Sen(wk, Nyear-1, NumYear, amssng, SenSlp,Prob,N1)
      SenSl_SPI_3m(ist,ilen) = SenSlp; Prob_SPI_3m(ist,ilen) = Prob
   Enddo
Enddo
End Subroutine Cal_Trend
!--------------------------------------------------------------
! Input: Mang chua SPI
! A. Tinh cho toan thoi ky:
! 1. Tan suat han: a) so thang han (SPI <= -1.5) hoac b) % so thang han tren tong so thang
! 2. Cuong do han: a) Tong cac SPI <= -1.5 hoac b) Trung binh SPI cac thang co SPI <= -1.5
! 3. Tong so dot han (Tong so dot co cac thang lien tiep SPI <= -1.5)
! 4. Do dai cac dot han 
! Ghi chu: 1-3: Moi o luoi/tram mot gia tri; 4: moi o luoi/tram la mot mang
! B. Tinh cho tung nam:
! 1. So thang han trong nam
! 2. Tong gia tri SPI <= -1.5 trong nam
! 3. ???
! Ghi chu: 
! 1) Tinh B1 roi tinh A1
! 2) A2b = A2a/A1a
! 3) Tinh A3 va A4 dong thoi
!--------------------------------------------------------------
Subroutine Cal_Drou_Characs
! Calculate all drought characteristics for all stations/gridpoints
Integer ist, ilen, iyr, imw
!Real :: Nguong
!Nguong = -1.5       ! Nguong han nang
!Do ist=1, Nst
!   Do ilen=1,NLen
!      print*,"Station=",Station(ist)
   !   print*,"Time Scale=",ilen
   ! For each year and each month first:
!      print*,"Calculating NDYrs(ist,1:NYear,ilen) ..."
      ! Determine Number of Drought months in each year
!      Call Cal_No_Drou_Mon_in_Yr (SPI(ist,1:NTime,ilen), NDYrs(ist,1:NYear,ilen))
!      print*,"Calculating SumSPIYrs(ist,1:NYear,ilen) ..."
!      ! Determine Magnitude (or severity) in each year (Total drought index in each year) 
!      Call Cal_Sum_Drou_IDX_in_Yr (SPI(ist,1:NTime,ilen), SumSPIYrs(ist,1:NYear,ilen))
!      print*,"Calculating DrEvent(ist,1:NYear,ilen) ..."

      ! Determine which month is Drought
!      Call Cal_Drou_Event_in_Whole_Period (SPI(ist,1:NTime,ilen), DrEvent(ist,1:Ntime,ilen))
!      print*,"Calculating Spell_Dura(ist,1:NYear,ilen) ..."

      ! Determine Drought Duration (k), frequency of each spell k-th, and number of spells
!      Call NSpell_and_Leng (DrEvent(ist,1:NTime,ilen), Ntime, Spell_Dura(ist,1:Ntime,ilen), Spell_Dura(ist,0,ilen))
           ! Spell_Dura(ist,0,ilen) la So DOT han cua tram ist ung voi qui mo ilen
           ! Spell_Dura(ist,k,ilen) (k=1,..,Spell_Dura(ist,0,ilen)) la DO KEO DAI cac dot  
!   Enddo
!Enddo
      print*,"Calculating DrEvent(ist,1:NYear,ilen) ..."
      ! Determine which month is Drought
      Call Cal_Drou_Event_in_Whole_Period ! (SPI(ist,1:NTime,ilen), DrEvent(ist,1:Ntime,ilen))

      print*,"Calculating Spell_Dura(ist,1:NYear,ilen) ..."
      ! Determine Drought Duration (k), frequency of each spell k-th, and number of spells
      Do ist=1, Nst
         Do ilen=1,NLen
            Call NSpell_and_Leng (DrEvent(ist,1:NTime,ilen), Ntime, Spell_Dura(ist,1:Ntime,ilen), Spell_Dura(ist,0,ilen))
         Enddo
      Enddo

      print*,"Calculating SumSPIYrs(ist,1:NYear,ilen) ..."
      ! Determine Magnitude (or severity) in each year (Total drought index in each year) 
      Call Cal_Sum_Drou_IDX_in_Yr !(SPI(ist,1:NTime,ilen), SumSPIYrs(ist,1:NYear,ilen))

      print*,"Calculating NDYrs(ist,1:NYear,ilen) ..."
      ! Determine Number of Drought months in each year
      Call Cal_No_Drou_Mon_in_Yr 
      ! For whole period:
      Print*,"Calculating ND_All_Yrs & SumSPI_All_Yrs."  
      Call Num_Drought_Mon_and_Sum_SPI
      
      Print*, "Calculating Geographycal Extention"
      Call Georaphical_Extent_1_Mon (Area_Ext1_2, 2) ! Area_Ext Thang 2
      Call Georaphical_Extent_1_Mon (Area_Ext1_3, 3) ! Area_Ext Thang 3
      Call Georaphical_Extent_1_Mon (Area_Ext1_4, 4) ! Area_Ext Thang 4

      Call Georaphical_Extent_3_Mon (Area_Ext3_234, 2)  ! 3 thang 2-3-4
      Call Magnitude_3_Mon (SPI_3m, 2)                  ! 3 thang 2-3-4

End Subroutine Cal_Drou_Characs
!--------------------------------------------------------------
Subroutine Num_Drought_Mon_and_Sum_SPI  ! Tong so thang han (ND_All_Yrs) va tong SPI cua cac thang han (SumSPI_All_Yrs) trong ca thoi ky
                                        ! (Chi tinh tu nam thu hai tro di)
Integer :: ist, ilen, it
Do ist=1, Nst
   Do ilen=1, NLen
      ND_All_Yrs(ist,ilen) = 0; SumSPI_All_Yrs(ist,ilen) = 0
      Do it=13, Ntime  ! 13 = Thang 1 nam thu 2
         if (DrEvent(ist,it,ilen) == 1) then
            ND_All_Yrs(ist,ilen) = ND_All_Yrs(ist,ilen) + DrEvent(ist,it,ilen)
            if (SPI(ist,it,ilen)<Nguong) SumSPI_All_Yrs(ist,ilen) = SumSPI_All_Yrs(ist,ilen) + SPI(ist,it,ilen)
         endif 
      Enddo
   Enddo
Enddo
End Subroutine Num_Drought_Mon_and_Sum_SPI
!--------------------------------------------------------------
Subroutine Georaphical_Extent_1_Mon (Area_Ext1, Mon)
! Tinh % tong so o luoi co han trong 1 thang so voi toan vung
Integer :: ist, ilen, it, Mon, iyr, N_VN, N_Dr
Real, Dimension(Nyear,NLen) :: Area_Ext1
Area_Ext1(:,:) = -99.
Do ilen=1, NLen
   iyr = 1  ! Bo nam dau tien, bat dau tu nam thu 2 tro di
   Do it=13, Ntime
      if (mod(it,12)==Mon) then  !.OR.mod(it,12)==0) then 
         iyr = iyr+1; N_VN = 0; N_Dr = 0
         Do ist=1, Nst
            if (VN_Inland (VN_Land(ist))) then
               N_VN = N_VN + 1
               if (DrEvent(ist,it,ilen)==1) then
                  N_Dr = N_Dr+1
               endif
             endif
         Enddo
         Area_Ext1(iyr,ilen) = Real(N_Dr)/Real(N_VN)*100
      endif
   Enddo
Enddo 
End Subroutine Georaphical_Extent_1_Mon
!--------------------------------------------------------------
Subroutine Georaphical_Extent_3_Mon (Area_Ext3, Mon)
! Tinh % tong so o luoi co han trong 3 thang lien tuc (Mon, Mon+1, Mon+2) so voi toan vung
!          >>>>>   Mon /= 12, Mon <=10 <<<<<
Integer :: ist, ilen, it, Mon, iyr, N_VN, N_Dr
Real, Dimension(Nyear,NLen) :: Area_Ext3
Area_Ext3(:,:) = -99.
Do ilen=1, NLen
   iyr = 1  ! Bo nam dau tien, bat dau tu nam thu 2 tro di
   Do it=13, Ntime-3   ! Phai co 3 thang lien tuc
      if (mod(it,12)==Mon) then
         iyr = iyr+1; N_VN = 0; N_Dr = 0
         Do ist=1, Nst
            if (VN_Inland (VN_Land(ist))) then
               N_VN = N_VN + 1
               if (DrEvent(ist,it,ilen)==1.AND.DrEvent(ist,it+1,ilen)==1.AND.DrEvent(ist,it+2,ilen)==1) then
                  N_Dr = N_Dr+1
               endif
             endif
         Enddo
         Area_Ext3(iyr,ilen) = Real(N_Dr)/Real(N_VN)*100
      endif
   Enddo
Enddo
End Subroutine Georaphical_Extent_3_Mon
!--------------------------------------------------------------
Subroutine Magnitude_3_Mon (SPI_3m, Mon)
! Tinh % tong chi so SPI 3 thang lien tuc (Mon, Mon+1, Mon+2) so voi toan vung
!          >>>>>   Mon /= 12, Mon <=10 <<<<<
Integer :: ist, ilen, it, Mon, iyr
Real :: Tmp
Real :: SPI_3m(Nst_Max,Nyear,NLen)
SPI_3m(:,:,:) = -99.
Do ist=1, Nst
   if (VN_Inland (VN_Land(ist))) then
      Do ilen=1, NLen
         if (.Not.Empty_Ar1(SPI(ist,1:Ntime,ilen), Ntime)) then
            iyr = 1  ! Bo nam dau tien, bat dau tu nam thu 2 tro di
            Do it=13, Ntime-3   ! Phai co 3 thang lien tuc
         !      Tmp = 0
               if (mod(it,12)==Mon) then
                  iyr=iyr+1
                  Tmp = SPI(ist,it,ilen)+SPI(ist,it+1,ilen)+SPI(ist,it+2,ilen)
                  SPI_3m(ist,iyr,ilen) = Tmp
               endif
         !      print*, SPI(ist,it,ilen)
         !      if (mod(it,12)==0.and.tmp/=0) then
         !         print*,VN_Land(ist)," ",iyr, " ", SPI_3m(ist,iyr,ilen)
         !         pause
         !      endif
            Enddo
         endif
      !   SPI_3m(ist,iyr,ilen) = Tmp
      Enddo
   endif
Enddo
End Subroutine Magnitude_3_Mon
!--------------------------------------------------------------
Logical Function VN_Inland (ist)
Real :: ist
if (ist==2.OR.ist==3.OR.ist==4.OR.ist==5.OR.ist==6.OR.ist==7.OR.ist==8) then
   VN_Inland = .True.
endif
End Function VN_Inland
!--------------------------------------------------------------
Character (Len=2) Function SubRegion (ist)
Real :: ist
if (ist==2) then
   SubRegion = "B1"
else if (ist==3) then
   SubRegion = "B2"
else if (ist==4) then
   SubRegion = "B3"
else if (ist==5) then
   SubRegion = "B4"
else if (ist==6) then
   SubRegion = "N1"
else if (ist==7) then
   SubRegion = "N2"
else if (ist==8) then
   SubRegion = "N3"
endif
End Function SubRegion
!--------------------------------------------------------------
Subroutine Cal_No_Drou_Mon_in_Yr      ! So thang han trong tung nam
! Tinh muc B1
! Input: Mang X(Ntime) la chi so SPI cua mot tram/o luoi; Nguong: Gia tri nguong xac dinh cuong do han
! Output: Mang NDYrs(Nyear) chua so thang han cua tung nam ung voi 1 qui mo thoi gian (1,3,6,12,... month)
Integer :: iyr, imw, i, ist, ilen
Real :: wk(Ntime)

!SPI(ist,1:NTime,ilen), NDYrs(ist,1:NYear,ilen)
Do ist=1, Nst
   Do ilen=1, Nlen
      wk(1:Ntime) = SPI(ist,1:NTime,ilen)
      if (Empty_Ar1 (wk, Ntime)) then
         NDYrs(ist,1:NYear,ilen) = amssng
      else
         i=0
         Do iyr=1, Nyear
            NDYrs(ist,iyr,ilen) = 0
            Do imw=1, NMon_Week
               i=i+1; if (wk(i)/=amssng .AND. wk(i)<Nguong) NDYrs(ist,iyr,ilen)=NDYrs(ist,iyr,ilen)+1
            Enddo
         Enddo
      endif
   Enddo
Enddo
End Subroutine Cal_No_Drou_Mon_in_Yr
!--------------------------------------------------------------
Subroutine Cal_Sum_Drou_IDX_in_Yr 
! Tinh muc B1
! Input: Mang X(Ntime) la chi so SPI cua mot tram/o luoi; Nguong: Gia tri nguong xac dinh cuong do han
! Output: Mang SumSPIYrs(Nyear) tong SPI<Nguong cua tung nam ung voi 1 qui mo thoi gian (1,3,6,12,... month)
! SPI(ist,1:NTime,ilen), SumSPIYrs(ist,1:NYear,ilen)
Integer :: iyr, imw, i, ist, ilen
Real :: wk(Ntime)

Do ist=1, Nst
   Do ilen=1, Nlen
      wk(1:Ntime) = SPI(ist,1:NTime,ilen)
      if (Empty_Ar1 (wk, Ntime)) then
         SumSPIYrs(ist,1:Nyear,ilen) = amssng
      else
         i=0
         Do iyr=1, Nyear
            SumSPIYrs(ist,iyr,ilen) = 0
            Do imw=1, NMon_Week
               i=i+1; if (wk(i)/=amssng) SumSPIYrs(ist,iyr,ilen)=SumSPIYrs(ist,iyr,ilen)+wk(i)  ! Tong SPI
              ! if (wk(i)/=amssng .AND. wk(i)<Nguong) SumSPIYrs(ist,iyr,ilen)=SumSPIYrs(ist,iyr,ilen)+wk(i)
            Enddo
         Enddo
      endif
   Enddo
Enddo
End Subroutine Cal_Sum_Drou_IDX_in_Yr
!--------------------------------------------------------------
Subroutine Cal_Drou_Event_in_Whole_Period !(X, Y)  ! Xac dinh cac thang co han (SPI < Nguong)
! Tinh muc B1
! Input: Mang X(Ntime) la chi so SPI cua mot tram/o luoi; Nguong: Gia tri nguong xac dinh cuong do han
! Output: Mang Y(Ntime) chua gia tri 0/1 cua SPI<Nguong ung voi 1 qui mo thoi gian (1,3,6,12,... month)
! (SPI(ist,1:NTime,ilen), DrEvent(ist,1:Ntime,ilen))
Integer :: iyr, imw, i, ist, ilen
Real :: wk(Ntime)

Do ist=1, Nst
   Do ilen=1, Nlen
      wk(1:Ntime) = SPI(ist,1:NTime,ilen)
      if (Empty_Ar1 (wk, Ntime)) then
         DrEvent(ist,1:Nyear,ilen) = amssng
      else
         Do i=1, Ntime
            if (wk(i) /= amssng .AND. wk(i) < Nguong) then
               DrEvent(ist,i,ilen) = 1
            else if (wk(i) /= amssng .AND. wk(i) >= Nguong) then
               DrEvent(ist,i,ilen) = 0
            else
               DrEvent(ist,i,ilen) = amssng
            endif
         Enddo
      endif
   Enddo
Enddo
End Subroutine Cal_Drou_Event_in_Whole_Period
!------------------------------------------------------------------------------
Subroutine NSpell_and_Leng (X1, N, M1, Imax1)   ! So dot han (Imax1) va do dai cac dot (M1(Imax1))
! Input: Array X(N) with values of 0, 1, or -99
! Output: Array M(Imax) with values of number of spell 
Integer :: N, i, k
!Integer :: X(NYear*NWeek), M(NYear*NWeek), Imax
Integer :: X(N), M(N), IMAX
Real :: X1(N), M1(N), IMAX1

!N = NTime ! NYear*NWeek
X(:) = X1(:)
!Print*,N

M(:) = 0; IMAX=0; I = 0; M1(:) = 0; IMAX1=0
DO
   I = I + 1  ! Loop until end of series data
   IF (I>N) EXIT
   IF (X(I)==1) THEN  ! Event occur at first time of this spell
      K = 0  ! Number of days event occur continuosly
      DO
         IF (X(I)==0.OR.X(I)==-99.OR.I==N) THEN   ! End of this spell
            IF (I==N.AND.X(I)==1) K=K+1   ! For the last spell
            IF (K/=0) M(K)=M(K)+1   ! Increase for this spell
            IMAX=IMAX+1             ! Increase number of spells
   !         PRINT*, K, M(K), IMAX
            EXIT
         END IF
         I = I + 1
         K = K + 1   ! How many days of this spell
      ENDDO
   ELSE
      CYCLE
   END IF
   IF (I>=N) EXIT
ENDDO
I=IMAX
!print*,"I, IMAX=",I, IMAX
if (I>0) then
   DO
      IF (M(I)/=0) EXIT    ! Determine exactly how many spells
      I=I-1
   ENDDO
   IMAX = I
! Gan lai cho so thuc:
   M1(:) = M(:); IMAX1 = IMAX
else
   M1(:) = 0; IMAX1 = 0
endif
!print*,"imax=",imax,"M=",M(1:imax)
End Subroutine NSpell_and_Leng
!--------------------------------End of this file--------------------------