! Chuong trinh tinh chi so han SPI da duoc kiem tra theo SL chuan
! Dong lenh goi Chuong trinh:
!     Call Cal_SPI
! Yeu cau bat buoc:
! Chuong trinh goi phai khai bao cac bien sau:
! 1) NYear: So nam co SL
! 2) Nmonth: So thang trong nam (=12)
! 3) Ntime = NYear*NMonth: Do dai chuoi (so thang)
! 4) NLen: So TimeScale se duoc tinh (vi du =4 de tinh cho han 1, 3, 6, 12 thang)
! 5) Len_SPI(NLen): Cac TimeScale (vi du Len_SPI(1)=1, Len_SPI(2)=3, Len_SPI(3)=6,...)
! 6) amssng: Missing value (vi du =-99)
! 7) Nst: So tram (hoac so diem luoi) se tinh
! 8) DatBL(1:Ntime,Nst): Mang chu SL mua thang
! 9) Rain(1:Ntime): Mang lam viec, chu SL mua thang
! 10) beta_spi(NMonth), gamm(NMonth), pzero(NMonth): Cac mang chua he so cac ham phan bo
! 11) SPI(Nst, Ntime, NLen): Chua ket qua tinh SPI cho Nst tram tu chuoi Ntime cho NLen TimeScale

!------------------------------------------------------------------------------
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/10) then ; Empty_Ar1 = .True. ; else ; Empty_Ar1 = .False. ; endif
End Function Empty_Ar1
!--------------------------------------------------------------
!------------------------------------------------------------------------------
Subroutine Cal_SPI
Integer :: ist, i
Loop_1: Do ist=1, Nst
   Rain(1:Ntime)=DatBL(1:Ntime,ist)
   if (Empty_Ar1(Rain, Ntime)) then
      spi(ist,:,:) = amssng
      Cycle Loop_1
   endif
 !  print "(12F7.1)",Rain(1:Ntime)
 !  pause
   Do i=1,NLen
      Call spigam (Len_SPI(i), Rain, beta_spi, gamm, pzero, spi(ist,:,i))
   Enddo
Enddo Loop_1
End Subroutine
!--------------------------------------------------------------
Subroutine spigam (nrun, pp, beta, gamm, pzero, idx)
! Input:
! nrun: Do dai dot han (1,3,6,12,... thang)
! pp: Chuoi luong mua thang xep theo thu tu thoi gian
! Output:
! idx: Mang 1 chieu chua chi so SPI (idx(:) se duoc gan cho spi(ist,:,i)

!Real :: pp(NYear*NMonth), idx(NYear*NMonth),tmparr(NYear+1), beta(NMonth), gamm(NMonth), pzero(NMonth)
Real :: pp(Ntime), idx(Ntime),tmparr(NYear+1), beta(Nmonth), gamm(NMonth), pzero(NMonth)
Integer :: nrun, j,i,n,im !, maxyrs
Real :: alpha
!   The first nrun-1 idx values will be missing.
  idx(1:nrun-1) = amssng
!  maxyrs = NYear
!
Loopj: Do j = nrun, maxyrs*Nmonth !NWeek !12
   idx(j) = 0.0
   Do i = 0, nrun-1
      if (pp(j - i) /= amssng) then
         idx(j) = idx(j) + pp(j - i)
      else
          idx(j) = amssng
          Cycle Loopj 
      endif
   Enddo
Enddo Loopj
Do i = 0, Nmonth - 1 !NWeek-1 ! 11
   n = 0
   Do j = nrun+i, maxyrs*Nmonth, Nmonth !NWeek, NWeek ! 12, 12
      if (idx(j) /= amssng) then
         n = n + 1; tmparr(n) = idx(j)
      endif
   Enddo
   im = mod (nrun+i-1, Nmonth) + 1 !NWeek) + 1 ! (nrun+i-1, 12) + 1
!     Here's where we do the fitting.
   Call gamfit (tmparr, n, alpha, beta(im), gamm(im), pzero(im))
Enddo
!     Replace precip. sums stored in index with SPI's
Do j = nrun, maxyrs*Nmonth !NWeek !12
   im = mod (j-1, Nmonth) + 1 !NWeek) + 1 ! (j-1,12) + 1
   if (idx(j) /= amssng) then
!     Get the probability
!     print*, im,beta(im), gamm(im), pzero(im)
      idx(j) = gamcdf(beta(im), gamm(im), pzero(im), idx(j))
!     Convert prob. to z value. 
      idx(j) = anvnrm(idx(j))
   endif
Enddo
End Subroutine 
!--------------------------------------------------------------
Real Function anvnrm (prob)
Real :: sgn,t,prob,c0, c1, c2,d1, d2, d3
data c0, c1, c2 /2.515517, 0.802853, 0.010328/
data d1, d2, d3 /1.432788, 0.189269, 0.001308/

      if (prob > 0.5) then
         sgn = 1.0
         prob = 1.0 - prob
      else
         sgn = -1.0
      endif

      if (prob < 0.0) then
         Print*, 'Error in anvnrm(). Prob. not in [0,1.0]'
         anvnrm = 0.0
         return
      endif

      if (prob == 0.0) then
         anvnrm = 1.0e37 * sgn
         return
      endif

      t = sqrt(alog (1.0 / (prob * prob)))
      anvnrm = (sgn * (t - ((((c2 * t) + c1) * t) + c0) / &
               ((((((d3 * t) + d2) * t) + d1) * t) + 1.0)))
      return
End Function
!--------------------------------------------------------------
Subroutine gamfit (datarr, n, alpha, beta, gamm, pzero)
Real :: datarr(*)
Real :: alpha, beta, gamm, pzero, sum,sumlog,av
Integer :: n, i,nact

      if (n <= 0) then
         Print*, 'Error in gamfit - empty data array'
         stop
      endif

      sum = 0.0
      sumlog = 0.0
      pzero = 0.0
      nact = 0

!     compute sums
Do i = 1, n
   if (datarr(i) > 0.0) then
      sum = sum + datarr(i)
      sumlog = sumlog + alog (datarr(i))
      nact = nact + 1
   else
      pzero = pzero + 1
   endif
Enddo
      pzero = pzero / n
      if(nact /= 0.0) av = sum / nact

!     Bogus data array but do something reasonable
      if (nact == 1) then
         alpha = 0.0
         gamm = 1.0
         beta = av
         return
      endif
!     They were all zeroes. 
      if (pzero == 1.0) then
         alpha = 0.0
         gamm = 1.0
         beta = av
         return
      endif

!     Use MLE
      alpha = alog (av) - sumlog / nact
      gamm = (1.0 + sqrt (1.0 + 4.0 * alpha / 3.0)) / (4.0 * alpha)
      beta = av / gamm

      return
End Subroutine
!--------------------------------------------------------------
Real function  gamcdf (beta, gamm, pzero, x)
Real :: beta, gamm, pzero, x
      if (x <= 0.0) then
         gamcdf = pzero
      else
         gamcdf = pzero + (1.0 - pzero) * gammap (gamm, x / beta)
      endif
End Function
!--------------------------------------------------------------
Real Function gaminv (beta, gamm, pzero, prob)
Real :: beta, gamm, pzero, prob, eps, thigh,tlow,phigh, t,p
Integer :: niter
data  eps /1.0e-7/

!     Check if prob < prob of zero
      if (prob <= pzero) then
         gaminv = 0.0
         return
      endif

!     Otherwise adjust prob
      prob = (prob - pzero) / (1.0 - pzero)

!     Make initial guess. Keep doubling estimate until prob is
!     bracketed.
      thigh = 2.0*eps
 Do 
    phigh = gamcdf (beta, gamm, pzero, thigh)
    if (phigh >= prob) exit
    thigh = thigh*2.0
 Enddo
      tlow = thigh / 2.0

!     Iterate to find root.
      niter = 0
 Do
    if ((thigh - tlow) <= eps) exit
    niter = niter + 1
    t = (tlow + thigh) / 2.0
    p = gamcdf (beta, gamm, pzero, t)
    if (p .lt. prob) then
       tlow = t
    else
       thigh = t
    endif
 Enddo
      gaminv = (tlow + thigh) / 2.0
End Function
!--------------------------------------------------------------
Real Function gamser (a, x)
Real :: a, x
!     Maximum number of iterations, and bound on error.
Real, parameter :: maxitr=100, eps=3.0e-7
Integer :: iwarn, n
Real :: gln,ap,del,sum
data iwarn /0/
 gln = gammln (a)
 if (x == 0.0) then
    gamser = 0.0
    return
 endif
 ap = a
 sum = 1.0 / a
 del = sum
 Do n = 1, maxitr
    ap = ap + 1.0
    del = del * (x / ap)
    sum = sum + del
    if (abs (del) < eps * abs (sum)) then !goto 20
       gamser =  sum * exp (-x + a * alog (x) - gln)
       return
    endif
 Enddo
 iwarn = iwarn + 1
 if (iwarn < 20) then
    Print*, 'gamser(',a,x,'): not converging.'
    Print*, 'Approximate value of ',sum,'  + /-',del,' used.'
 endif
 gamser =  sum * exp (-x + a * alog (x) - gln)
End Function
!--------------------------------------------------------------
Real Function gammcf (a, x)
Real, parameter :: eps=3.0e-7
Integer, parameter :: maxitr=100
Real :: a, x, a0,a1,b0,b1,fac,gln,gold,anf,an,ana,g
Integer :: nwarn, n
Data nwarn / 0 /, g / 0.0 /

 gln = gammln (a)
 gold = 0.0
 a0 = 1.0
 a1 = x
 b0 = 0.0
 b1 = 1.0
 fac = 1.0
 Do n = 1, maxitr
    an = n
    ana = an - a
    a0 = (a1 + a0 * ana) * fac
    b0 = (b1 + b0 * ana) * fac
    anf = an * fac
    a1 = x * a0 + anf * a1
    b1 = x * b0 + anf * b1
    if (a1 /= 0.0) then
       fac = 1.0 / a1
       g = b1 * fac
       if (abs((g - gold) / g) .lt. eps) then !goto 20
          gammcf =  g * exp (-x + a * alog (x) - gln)
          return
       endif
       gold = g
    endif
 Enddo
 nwarn = nwarn + 1
 if (nwarn < 20) then
    Print*, 'gammcf(',a,x,'): not converging.'
    Print*, 'Inaccurate value of ', g, ' +/- ', abs(g - gold), ' used.'
 endif
 gammcf =  g * exp (-x + a * alog (x) - gln)
End Function
!--------------------------------------------------------------
Real Function gammap (a, x)
Real :: a, x
      if (x < a + 1.0) then
         gammap = gamser (a, x)
      else
         gammap = 1.0 - gammcf (a, x)
      endif
End Function
!--------------------------------------------------------------
Real Function gammaq (a, x)
Real :: a, x
      if (x < a + 1.0) then
         gammaq = 1.0 - gamser (a, x)
      else
         gammaq = gammcf (a, x)
      endif
End Function
!--------------------------------------------------------------
Real Function gammln(xx)
Real :: x, xx, tmp, ser
Real :: cof(6)
data cof /76.18009173, -86.50532033, 24.01409822, -1.231739516, 0.120858003e-2, -0.536382e-5/
      x = xx - 1.0
      tmp = x + 5.5
      tmp = tmp - (x+0.5) * alog (tmp)
      ser = 1.0
      Do j = 1, 5
         x = x + 1.0
         ser = ser + cof(j) / x
      Enddo
      gammln = -tmp + alog (2.50662827465 * ser)
End Function
!--------------------------------------------------------------