PROGRAM SPEI
  implicit none
  integer, parameter :: NUMDATOSMAX = 5000, NUMRESULTMAX = 5000, NUMSEASONSMAX = 12
  character (len = 100) :: pathOrigen, pathDestino
  character (len =  20) :: estacion, latitud
  character (len =   2) :: c_acumulated
  integer :: io
  integer :: anio, mes, seasonality, acumulated, numRegistros, acumRegistros, indice, jndice
  real :: lat
  real, dimension (0:NUMDATOSMAX-1) :: rainSeries,    tempSeries, etpSeries, &
                                       balanceSeries, acumSeries, speiSeries    

  !! ========================================================================================
  !! Initialize variables

  anio         = 0 ; mes           = 0 ; seasonality = 0 ; acumulated = 0
  numRegistros = 0 ; acumRegistros = 0 ; indice      = 0 ; jndice     = 0  

  lat          = 0.
  rainSeries    (:) = 0. ; tempSeries (:) = 0. ; etpSeries  (:) = 0.
  balanceSeries (:) = 0. ; acumSeries (:) = 0. ; speiSeries (:) = 0.

  !! ----------------------------------------------------------------------------------------
  !! Read in-line arguments; Ex: ./spei.exe nlength input.txt output.txt

  Call getarg (1,c_acumulated)  ; read (c_acumulated,*) acumulated
  Call getarg (2,pathOrigen)
  Call getarg (3,pathDestino)

  !! ----------------------------------------------------------------------------------------
  !! Read data

  print*, "- Reading file: ", pathOrigen
  open (1, file = pathOrigen)
    read (1,*) estacion         ! Ten tram
    read (1,*) lat              ! Latitude
    read (1,*) anio, mes        ! Nam bat dau, Thang bat dau (cua chuoi thoi gian)
    read (1,*) seasonality      ! Chu ky mua (1 nam - 12 thang)
    if (seasonality > NUMSEASONSMAX) stop !! Maximum is 12 months (NO FIX)
    
    indice = 0                            !! Array index start from 0
    do 
      if (indice > NUMDATOSMAX) stop      !! Maximum is 5000 line (can be fixed above)
      
      read (1,*, iostat = io) rainSeries(indice), tempSeries(indice)  ! Doc SL Preci, Temp
      
      if (io .ne. 0) exit
      indice = indice + 1
    enddo
    numRegistros = indice
    if ( tempSeries(1) == 0.) numRegistros = numRegistros - 1   !! ???? 
    print*, "NumRegistros (Do dai chuoi - So thang)= ",numRegistros
  close (1)

  !! Print check
  print*, "- Station     : ", estacion
  print*, "- Latitude    : ", lat
  print*, "- Initial date: ", mes, anio
  print*, "- Seasonality : ", seasonality
  print*, "- Register    : ", numRegistros
  print*, "- Period      : ", acumulated
  print*, "---------------------------------------"

  !! ----------------------------------------------------------------------------------------
  !! Compute the climatic balance: precipitation minus potential evapotranspiration

!  if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) ) then   ! Kiem tra chuoi Nhiet do
!    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
!    do indice = 0, numRegistros - 1
!      balanceSeries(indice) = rainSeries(indice)    ! Bo qua boc hoi neu khong co nhiet do
!    enddo
!  endif

  !! ----------------------------------------------------------------------------------------
  !! Compute the cumulative series  ! 

  !anio = anio + nint(real((acumulated - 1) / 12))

  ! acumulated: Quy mo thoi gian cua han (Drought Time-scale)

!  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

  !! ----------------------------------------------------------------------------------------
  !! Compute the SPEI series

  Call Calculate_Water_Balance (lat, tempSeries, rainSeries, balanceSeries, numRegistros)
  Call Calculate_Accumulated_Series (acumulated, balanceSeries, numRegistros, acumSeries, acumRegistros)

  Call spei(acumSeries, acumRegistros, seasonality, speiSeries)

  !! ----------------------------------------------------------------------------------------
  !! Write result to file

  print*, "- Writing to file: ", pathDestino
  open (2, file = pathDestino )
    write (2,*) estacion
    write (2,*) lat
    write (2,*) anio,mes
    write (2,*) seasonality
    print*,"numRegistros,acumRegistros,acumulated=",numRegistros,acumRegistros,acumulated
    do indice = 0, acumRegistros - 1
      write (2,*) speiSeries(indice)
    enddo
  close (2)

!! ==========================================================================================
CONTAINS
!---------------------------------------------------------------------------------!
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, numRegistros, acumRegistros
Real, Dimension (0:NUMDATOSMAX-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 "(A,10I5)","numRegistros,acumulated,indice,jndice=", numRegistros,acumulated,indice,jndice
      !! print*,"CHECKING:", acumSeries(indice - acumulated + 1) !! NXTTNX
    enddo
!    pause
  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:NUMDATOSMAX-1) :: tempSeries, rainSeries, etpSeries, balanceSeries
Real :: lat

  if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) ) then   ! Kiem tra chuoi Nhiet do
    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
    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
!---------------------------------------------------------------------------------!
Subroutine spei(dataSeries, n, seasons, speiSeries)
  integer :: n, seasons
  real, dimension (0:n-1) :: dataSeries, speiSeries
  integer :: i, j, k, nSeason
  real :: seasonSeries (0:NUMDATOSMAX-1), beta (0:2), logLogisticParams(0:NUMSEASONSMAX-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"
END