!------------------------------------------------------------------------------
Subroutine Cal_SPEI
Integer :: ist, i, acumRegistros, TimeScale
Real, Dimension (0:Ntime-1) :: Rain, T2m, Rain0, T2m0
Real, Dimension (0:Ntime-1) :: balanceSeries, acumSeries, speiSeries

Loop_1: Do ist=1, Nst
   ! Set to Zero before calculating:
   Rain0 (:) = 0.; T2m0 (:) = 0.; spi(ist,:,:) = 0.
!   balanceSeries (:) = 0. ; acumSeries (:) = 0.; speiSeries (:) = 0.
   ! Now calcutale:
   Rain0(0:Ntime-1)=DatBL(1:Ntime,ist); T2m0(0:Ntime-1)=DatTMP(1:Ntime,ist)
   if (Empty_Ar1(Rain, Ntime-1).OR.Empty_Ar1(T2m, Ntime-1)) then
      spi(ist,:,:) = amssng
      Cycle Loop_1
   endif
!   print "(12F7.1)",Rain(0:Ntime-1)
!   print "(12F7.1)",T2m(0:Ntime-1)
!   pause
   Do i=1,NLen
      balanceSeries (:) = 0. ; acumSeries (:) = 0.; speiSeries (:) = 0.;
      Rain(:)=Rain0(:); T2m(:)=T2m0(:)   ! Luu lai cho lan tinh sau
      TimeScale = Len_SPI(i)
!      Call Calculate_Water_Balance (lat(ist), T2m, Rain, balanceSeries, Ntime)
!      Call Calculate_Accumulated_Series (TimeScale, balanceSeries, Ntime, acumSeries, acumRegistros)
!      Call spei(acumSeries, acumRegistros, NMonth, speiSeries)
!      Call spei(acumSeries, acumRegistros, TimeScale, speiSeries)
      Call Calc_acumSeries (T2m, Rain, Lat(ist), Ntime, acumRegistros, TimeScale, speiSeries)
      spi(ist,TimeScale:Ntime,i) = speiSeries(0:acumRegistros-1)
!      if (i==1.And.ist>1) then; print "(12F8.3)",spi(ist,1:Ntime, i)-spi(ist-1,1:Ntime, i); pause; endif
   Enddo
 !  print "(F12.3)", lat(ist)
 !  print "(12F8.1)", T2m(0:Ntime-1); print "(12F8.1)", Rain(0:Ntime-1); print "(12F10.3)", spi(ist,1:Ntime, NLen); pause
Enddo Loop_1
End Subroutine
!---------------------------------------------------------------------------------!
Subroutine Calc_acumSeries (tempSeries, rainSeries, lat, numRegistros, acumRegistros, acumulated, speiSeries) 
  integer :: acumulated, numRegistros, acumRegistros, indice, jndice
  real :: lat
  real, dimension (0:Ntime-1) :: rainSeries, tempSeries, etpSeries, balanceSeries, acumSeries, speiSeries

  !! Compute the climatic balance: precipitation minus potential evapotranspiration
  if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) ) then
    Call thornthwaite (tempSeries, numRegistros, lat, etpSeries)
    do indice = 0, numRegistros - 1
      balanceSeries(indice) = rainSeries(indice) - etpSeries(indice)
    enddo
  else
    do indice = 0, numRegistros - 1
      balanceSeries(indice) = rainSeries(indice)
    enddo
  endif
!  print "(12F8.3)", tempSeries(0:Ntime-1), balanceSeries(0:Ntime-1); pause;
  !----------------------------------------------------------------------------------------
  !! Compute the cumulative series
  acumRegistros = numRegistros - acumulated + 1; acumSeries(:) = 0.
  do indice = acumulated - 1, numRegistros -1
    do jndice = 0, acumulated - 1
      acumSeries(indice - acumulated + 1) = &
                  acumSeries(indice - acumulated + 1) + balanceSeries (indice - jndice)
    enddo
  enddo
  !----------------------------------------------------------------------------------------
  !! Compute the SPEI series
  Call spei(acumSeries, acumRegistros, NMonth, speiSeries)
!  print "(12F12.3)", speiSeries(0:50); pause;
  
End Subroutine Calc_acumSeries
!---------------------------------------------------------------------------------!
Subroutine Calculate_Accumulated_Series (acumulated, balanceSeries, numRegistros, acumSeries, acumRegistros)
  ! acumulated: Quy mo thoi gian cua han (Drought Time-scale)
  ! Tinh chuoi tong tich luy CAN BANG NUOC (balanceSeries) trong thoi doan quy mo han (acumulated)
! Input: acumulated, balanceSeries(0:numRegistros-1)
! Output: acumSeries(0:numRegistros-1), acumRegistros

Integer :: acumulated, indice, jndice, numRegistros, acumRegistros
Real, Dimension (0:Ntime-1) :: balanceSeries, acumSeries

!  anio = anio + (acumulated - 1) / 12
!  mes  = mes  + acumulated - 1
!  print*, "anio, mes=", anio, mes; pause

!  do while (mes > 12)
!    mes = mes - 12
!  enddo

  acumRegistros = numRegistros - acumulated + 1

!  print*, "mes, acumRegistros=",  mes, acumRegistros; pause

  do indice = acumulated - 1, numRegistros -1
    do jndice = 0, acumulated - 1
      acumSeries(indice - acumulated + 1) = &
                  acumSeries(indice - acumulated + 1) + balanceSeries (indice - jndice)
      !! print*,"CHECKING:", acumSeries(indice - acumulated + 1) !! NXTTNX
    enddo
  enddo
End Subroutine Calculate_Accumulated_Series
!---------------------------------------------------------------------------------!
Subroutine Calculate_Water_Balance (lat, tempSeries, rainSeries, balanceSeries, numRegistros) 
!! Compute the climatic balance: precipitation minus potential evapotranspiration
! Input: Lat, tempSeries(0:numRegistros-1)
! Output: balanceSeries(0:numRegistros-1)
Integer :: indice, numRegistros
Real, Dimension (0:Ntime-1) :: tempSeries, rainSeries, etpSeries, balanceSeries
Real :: lat
  etpSeries (:) = 0
  if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) ) then   ! Kiem tra chuoi Nhiet do
!   print*,"if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) )"
!   print*, tempSeries(1), tempSeries(2)
    Call thornthwaite (tempSeries, numRegistros, lat, etpSeries)      ! Tinh boc hoi (etpSeries)
    do indice = 0, numRegistros - 1
      balanceSeries(indice) = rainSeries(indice) - etpSeries(indice)
      !! print*,"CHECKING", rainSeries(indice), etpSeries(indice), balanceSeries(indice)
    enddo
  else
    print*,"Bo qua boc hoi neu khong co nhiet do"
    do indice = 0, numRegistros - 1
      balanceSeries(indice) = rainSeries(indice)    ! Bo qua boc hoi neu khong co nhiet do
    enddo
  endif
End Subroutine Calculate_Water_Balance
!------------------------------------------------------------------------------
Logical Function Empty_Ar1 (X, N)
! Tra ve .True. neu X(N) khuyet toan bo SL
Integer :: N, i, k
Real :: X(0:N)
k = 0 ! Dem so phan tu khuyet cua X
Do i=0, N
   if (X(i) == amssng) k = k + 1
Enddo
if (k > N/10) then ; Empty_Ar1 = .True. ; else ; Empty_Ar1 = .False. ; endif
End Function Empty_Ar1
!---------------------------------------------------------------------------------!
Subroutine spei(dataSeries, n, seasons, speiSeries)
  integer :: n, seasons
!  real, dimension (0:n-1) :: dataSeries, speiSeries
  real, dimension (0:Ntime-1) :: dataSeries, speiSeries
  integer :: i, j, k, nSeason
  real :: seasonSeries (0:Ntime-1), beta (0:2), logLogisticParams(0:Ntime-1,0:2)

  do j = 1, seasons
    k = 0
    do i = j - 1, n - 1, seasons
      seasonSeries (k) = dataSeries (i)
      k = k + 1
    enddo
    nSeason = k
    Call upward (seasonSeries, nSeason)

    !! Compute probability weighted moments
    Call pwm (seasonSeries, nSeason, beta, 0., 0., 0)
    !!print*,"beta",beta(0), beta(1), beta(2)

    !! Fit a Log Logistic probability function
    Call logLogisticFit (beta, logLogisticParams(j,:))
    !! print*,logLogisticParams(j,0), logLogisticParams(j,1), logLogisticParams(j,2)

    !! Calculate the standardized values
    do i = j - 1, n - 1, seasons
      speiSeries (i) = logLogisticCDF (dataSeries (i), logLogisticParams (j,:))
      !print*,standardGaussianInvCDF (speiSeries(i))
      speiSeries (i) = -standardGaussianInvCDF (speiSeries(i))
      !print*,"CHECK: ",speiSeries (i)
    enddo
  enddo
End Subroutine spei
!! --------------------------
include "auxiliary.F90"
include "thornthwaite.F90"
include "lmoments.F90"
include "pdfs.F90"