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
    read (1,*) lat
    read (1,*) anio, mes
    read (1,*) seasonality
    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)
      
      if (io .ne. 0) exit
      indice = indice + 1
    enddo
    numRegistros = indice
    if ( tempSeries(1) == 0.) numRegistros = numRegistros - 1   !! ???? 
  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
    Call thornthwaite (tempSeries, numRegistros, lat, 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)
    enddo
  endif

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

  !anio = anio + nint(real((acumulated - 1) / 12))
  anio = anio + (acumulated - 1) / 12
  mes  = mes  + acumulated - 1
  do while (mes > 12)
    mes = mes - 12
  enddo
  acumRegistros = numRegistros - acumulated + 1
  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 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
    do indice = 0, acumRegistros - 1
      write (2,*) speiSeries(indice)
    enddo
  close (2)

!! ==========================================================================================
CONTAINS
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"
!auxiliary.F90  lmoments.F90  pdfs.F90  spei.F90  thornthwaite.F90
END