! Chuong trinh tinh chi so han SPI da duoc kiem tra theo SL chuan
! Dong lenh goi Chuong trinh:
!     Call Cal_SPEI
! Yeu cau bat buoc:
! Chuong trinh goi phai khai bao cac bien sau:
! 1) Nmonth: So thang trong nam (=12)
! 2) Ntime: Do dai chuoi (so thang) = NYear*NMonth
! 3) NLen: So TimeScale se duoc tinh (vi du =4 de tinh cho han 1, 3, 6, 12 thang)
! 4) Len_SPI(NLen): Cac TimeScale (vi du Len_SPI(1)=1, Len_SPI(2)=3, Len_SPI(3)=6,...)
! 5) amssng: Missing value (vi du =-99)
! 6) Nst: So tram (hoac so diem luoi) se tinh

! 7) DatPr(1:Ntime,Nst), DatTx(1:Ntime,Nst), DatTm(1:Ntime,Nst), DatT2m(1:Ntime,Nst) : Mang chu SL mua thang cua Nst tram

! 8) SPI(Nst, Ntime, NLen): Chua ket qua tinh SPI cho Nst tram tu chuoi Ntime cho NLen TimeScale
! 9) SPEI_Pars(Nst,0:NMonth,0:2,NLen)
!10) LogParams(0:NMonth,0:2)
!------------------------------------------------------------------------------

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

Loop_1: Do ist=1, Nst
   ! Set to Zero before calculating:
   Rain0 (:) = 0.; T2m0 (:) = 0.; 
   ! Now calcutale:
   !!! Rain0(0:Ntime-1)=DatBL(1:Ntime,ist); T2m0(0:Ntime-1)=DatTMP(1:Ntime,ist)
   !! Rain0(0:Ntime-1)=DatPr(1:Ntime,ist); Tx0(0:Ntime-1)=DatTx(1:Ntime,ist); Tm0(0:Ntime-1)=DatTm(1:Ntime,ist)
   Rain0(0:Ntime-1)=DatPr(1:Ntime,ist); T2m0(0:Ntime-1)=DatT2m(1:Ntime,ist)
   if (Empty_Arr1(Rain0, Ntime-1).OR.Empty_Arr1(T2m0, Ntime-1)) then
   !! if (Empty_Arr1(Rain0, Ntime-1).OR.Empty_Arr1(Tx0, Ntime-1).OR.Empty_Arr1(Tm0, Ntime-1)) then
      spi(ist,:,:) = amssng
      Cycle Loop_1
   endif
   print*,"SPEI - Calculating station: ", ist, Ntime
   !if (ist>=31950) print*,Rain, T2m
   Do i=1,NLen
      balanceSeries (:) = 0. ; acumSeries (:) = 0.; 
      Rain(:)=Rain0(:); T2m(:)=T2m0(:)   ! Luu lai cho lan tinh sau
      !Tx(:)=Tx0(:); Tm(:)=Tm0(:)
      TimeScale = Len_SPI(i)
      !if (ist>=32460) then
      !   print*,Rain
      !   print*,T2m
      !endif
      !print*,TimeScale
      Call Calc_SPEI_thornthwaite (T2m, Rain, Lat(ist), Ntime, acumRegistros, TimeScale, speiSeries)
      !!Call Calc_SPEI_Idx (Tx, Tm, Rain, Lat(ist), Ntime, acumRegistros, TimeScale, speiSeries)
      spi(ist,TimeScale:Ntime,i) = speiSeries(0:acumRegistros-1)
      !if (ist>=31950) print*, SPI

      SPEI_Pars(ist,0:NMonth,0:2,i) = LogParams(0:NMonth,0:2)

   Enddo
Enddo Loop_1
End Subroutine Cal_SPEI_thornthwaite
!---------------------------------------------------------------------------------!
Subroutine Cal_SPEI_Hargreaves
Integer :: ist, i, acumRegistros, TimeScale
!!! Real, Dimension (0:Ntime-1) :: Rain, T2m, Rain0, T2m0
Real, Dimension (0:Ntime-1) :: Rain, Tx, Tm, Rain0, Tx0, Tm0
Real, Dimension (0:Ntime-1) :: balanceSeries, acumSeries, speiSeries

Loop_1: Do ist=1, Nst
   ! Set to Zero before calculating:
   Rain0 (:) = 0.; Tx(:)=0.; Tm(:)=0; !T2m0 (:) = 0.;
   ! Now calcutale:
   !!! Rain0(0:Ntime-1)=DatBL(1:Ntime,ist); T2m0(0:Ntime-1)=DatTMP(1:Ntime,ist)
   Rain0(0:Ntime-1)=DatPr(1:Ntime,ist); Tx0(0:Ntime-1)=DatTx(1:Ntime,ist); Tm0(0:Ntime-1)=DatTm(1:Ntime,ist)
   !!! if (Empty_Arr1(Rain0, Ntime-1).OR.Empty_Arr1(T2m0, Ntime-1)) then
   if (Empty_Arr1(Rain0, Ntime-1).OR.Empty_Arr1(Tx0, Ntime-1).OR.Empty_Arr1(Tm0, Ntime-1)) then
      spi(ist,:,:) = amssng
      Cycle Loop_1
   endif
   print*,"Calculating station: ", ist
   !if (ist>=31950) print*,Rain, T2m
   Do i=1,NLen
      balanceSeries (:) = 0. ; acumSeries (:) = 0.;
      Rain(:)=Rain0(:); ! T2m(:)=T2m0(:)   ! Luu lai cho lan tinh sau
      Tx(:)=Tx0(:); Tm(:)=Tm0(:)
      TimeScale = Len_SPI(i)
      !if (ist>=32460) then
      !   print*,Rain
      !   print*,T2m
      !endif
      !print*,TimeScale
      !!! Call Calc_SPEI_Idx (T2m, Rain, Lat(ist), Ntime, acumRegistros, TimeScale, speiSeries)
      Call Calc_SPEI_Hargreaves (Tx, Tm, Rain, Lat(ist), Ntime, acumRegistros, TimeScale, speiSeries)
      spi(ist,TimeScale:Ntime,i) = speiSeries(0:acumRegistros-1)
      !if (ist>=31950) print*, SPI

      SPEI_Pars(ist,0:NMonth,0:2,i) = LogParams(0:NMonth,0:2)

   Enddo
Enddo Loop_1
End Subroutine Cal_SPEI_Hargreaves
!---------------------------------------------------------------------------------!
!!! Subroutine Calc_SPEI_Idx (tempSeries, rainSeries, lat, numRegistros, acumRegistros, acumulated, speiSeries) 
Subroutine Calc_SPEI_Hargreaves (Tx, Tm, 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
  real, dimension (0:Ntime-1) :: rainSeries, Tx, Tm, 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)
  if (Tx(1)/=0.AND.Tx(2)/=0.AND.Tm(1)/=0.AND.Tm(2)/=0) then
    Call Hargreaves (Tx, Tm, 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
  !----------------------------------------------------------------------------------------
  !! 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)
  
End Subroutine Calc_SPEI_Hargreaves
!---------------------------------------------------------------------------------!
Subroutine Calc_SPEI_thornthwaite (tempSeries, rainSeries, lat, numRegistros, acumRegistros, acumulated, speiSeries) 
  ! numRegistros = Ntime
  ! acumulated   = Nlen (Timescale)
  ! acumRegistros= Ntime - Nlen (Do dai chuoi tong truot)
  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
  !----------------------------------------------------------------------------------------
  ! 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)

End Subroutine Calc_SPEI_thornthwaite
!------------------------------------------------------------------------------
Logical Function Empty_Arr1 (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
     !if (X(i)<=0) k = k + 1
  Enddo
  if (k > N/10) then ; Empty_Arr1 = .True. ; else ; Empty_Arr1 = .False. ; endif
  !if (k > 0) then ; Empty_Arr1 = .True. ; else ; Empty_Arr1 = .False. ; endif
End Function Empty_Arr1
!---------------------------------------------------------------------------------!
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:NMonth,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
  LogParams(0:NMonth,0:2) = logLogisticParams(0:NMonth,0:2)
End Subroutine spei
!! --------------------------
!include "auxiliary.F90"
!include "thornthwaite.F90"
!include "lmoments.F90"
!include "pdfs.F90"
!=============================================================
! auxiliary.F90:
!! Returns the natural logarithm of the gamma function, ln[gamma(xx)] for xx > 0
real Function gammaLn (xx)
  real :: xx
  real :: x, y, tmp, ser
  real :: cof (0:5) = (/76.18009172947146, -86.50532032941677, 24.01409824083091,&
                        -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5/)
  integer :: j

  y = xx ; x = xx
  tmp = x + 5.5
  tmp = tmp - ((x + 0.5) * log(tmp))
  ser = 1.000000000190015
  do j = 0, 5
    y = y + 1
    ser = ser + cof(j) / y
  enddo

  gammaLn = - tmp + log (2.5066282746310005*ser / x)

End Function gammaLn
!------------------------------------------------------------------------------
!! Sorts a given data series from the lowest to the highest value
Subroutine upward (series, n)
  integer :: n
  real :: series (0:Ntime-1)  !series (0:n-1)
  integer :: i, j
  real :: temp

  do i = 0, n - 2
    do j = 0, n - 2 - i
      if ( series (j + 1) .lt. series (j) ) then
        temp = series (j)
        series (j) = series (j + 1)
        series (j + 1) = temp
      endif
    enddo
  enddo

!  do i = 0, n-1
!    print*, series (i)
!  enddo
End Subroutine upward
!-------------------------------------------------------------------
! lmoments.F90:
!! Calculates the first three probability weighted moments of a sample,
!! using either the unbiased estimator (when A=B=0)
!! or a plotting position formula (when A<=B<=0).
!! This are alpha PWMs, following Rao & Hamed 2000, eqs. 3.1.4 and 3.1.6
Subroutine pwm (series, n, pwms, A, B, isBeta)
  integer :: n, isBeta
  real :: series (0:n-1), pwms(0:2), A, B
  integer :: i
  real :: acum(0:2), F

  acum (:) = 0.
  if ( (A .eq. 0.) .and. (B .eq. 0.) ) then     !! use unbiased estimator
    do i = 1, n
      acum(0) = acum (0) + series (i-1)
      if (isBeta .eq. 0) then                   !! compute alpha PWMs
        acum(1) = acum(1) + series (i - 1) * (n - i)
        acum(2) = acum(2) + series (i - 1) * (n - i) * (n - i - 1)
      endif

      if (isBeta .eq. 1) then                   !! compute beta PWMs
        acum(1) = acum(1) + series (i - 1) * (i - 1)
        acum(2) = acum(2) + series (i - 1) * (i - 1) * (i - 2)
      endif
    enddo
  else  !! use plotting-position (biased) estimator
    do i = 1, n
      acum(0) = acum (0) + series (i-1)
      F = (real(i) + A) / (real(n) + B)
      if (isBeta .eq. 0) then                   !! compute alpha PWMs
        acum (1) = acum (1) + series(i - 1) * (1 - F)
        acum (2) = acum (2) + series(i - 1) * (1 - F) * (1 - F)
      endif
      if (isBeta == 1) then                     !! compute beta PWMs
        acum (1) = acum (1) + series (i - 1) * F
        acum (2) = acum (2) + series (i - 1) * F * F
      endif
    enddo
  endif

  pwms (0) = acum (0) / n
  pwms (1) = acum (1) / n / (n - 1)
  pwms (2) = acum (2) / n / ((n-1)*(n-2))

  return
End Subroutine pwm

!-------------------------------------------------------------------
! pdfs.F90: Estimates the parameters of a Gamma distribution functions
Subroutine logLogisticFit (beta, logLogisticParams)
  real, dimension (0:2) :: beta, logLogisticParams
  real :: g1, g2

  !! estimate gamma parameter
  logLogisticParams (2) = (2*beta (1) - beta (0)) / (6*beta(1)-beta(0)-6*beta(2))
  g1 = exp(gammaLn(1+1/logLogisticParams(2)))
  g2 = exp(gammaLn(1-1/logLogisticParams(2)))
  !! estimate alpha parameter
  logLogisticParams(1) = (beta(0)-2*beta(1))*logLogisticParams(2) / (g1*g2)
  !! estimate beta parameter
  logLogisticParams(0) = beta(0) - logLogisticParams(1)*g1*g2

  !!print*,"beta",beta(0), beta(1), beta(2)
  !!print*,"LOG ",logLogisticParams(0), logLogisticParams(1),
  !logLogisticParams(2)
  return
End Subroutine logLogisticFit
!------------------------------------------------------------------------------
!! Gives the cumulative distribution function of 'value', following a LogLogistic distribution
real Function logLogisticCDF (val, params)
  real :: val, params (0:2)
  logLogisticCDF = 1 / (1 + (params(1)/(val-params(0)))**params(2))
  return
end Function logLogisticCDF
!------------------------------------------------------------------------------
!! Gives the inverse cumulative distribution function following a standard Gaussian
!! distribution. I.e., the function gives the standardized value (Z) corresponding to a
!! given probability (0<prob<1) following a standard Gaussian distribution (mean=0 and sd=1).

real Function standardGaussianInvCDF (prob)
  real :: prob
  real :: C(0:2) = (/ 2.515517, 0.802853, 0.010328 /)
  real :: d(0:3) = (/ 0.000000, 1.432788, 0.189269,0.001308/)
  real :: W, WW, WWW, resul

  if (prob .le. 0.5) then
    W = sqrt(-2*log(prob))
  else
    W = sqrt(-2*log(1-prob))
  endif
  WW = W*W
  WWW = WW*W
  resul = W - (C(0) + C(1)*W + C(2)*WW) / (1 + d(1)*W + d(2)*WW + d(3)*WWW)

  if (prob .gt. 0.5) resul = -resul
  !print*, resul

  standardGaussianInvCDF = resul
  return
end Function standardGaussianInvCDF

!-------------------------------------------------------------------
! thornthwaite.F90:
Subroutine thornthwaite (tempSeries, n, lat, etpSeries)
  !! Calculates the potential evapotranspiration from a series of monthly
  !! temperature data, following the method of Thornthwaite.
  integer :: n
  real :: tempSeries (0:n-1), lat, etpSeries (0:n-1)

  real :: T, J, J2, J3, c, N_up, omega, K_up(13), tanLatMonth
  integer :: month, i, k

  real, parameter ::  pi = 3.14159265358979
  !! Tangens of the average solar declination angle for each month of the year
  real, parameter ::  tanDelta (0:11) = (/-0.37012566, -0.23853358, -0.04679872, 0.16321764, &
                                           0.32930908,  0.40677729,  0.37477410, 0.23906300, &
                                           0.04044485, -0.16905776, -0.33306377, -0.40743608/)
  real, parameter ::  days (0:11) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

  real :: tanLat
  tanLat = tan(0.0174532925*lat)

  !print*,"tanLat:", tanLat
  !!
  !-----------------------------------------------------------------------------------------
  !! Compute J: annual temperature efficiency index
  J = 0.
  do month = 1, 12
    T = 0. ; k = 0
    do i = month - 1, n - 1, 12
      T = T + tempSeries (i)
      k = k + 1
    enddo
    !! T: monthly average temperature
    T = T / real(k)
    !! print*,"CHECK:", J !! NXTTNX
    if (T > 0.) J = J + (T / 5.)**1.514
  enddo

  !! Compute c exponent
  J2 = J*J
  J3 = J2*J
  c = 0.000000675*J3 - 0.0000771*J2 + 0.01792*J + 0.49239

  !! print*,"CHECK:", c !! NXTTNX
  !! Compute K: monthly correction factor, depending on latitude
  do month = 1, 12
    tanLatMonth = tanLat*tanDelta(month-1)
    if ( (tanLatMonth .gt. -1.) .and. (tanLatMonth .lt. 1.) ) omega = acos(-tanLatMonth)
    if ( tanLatMonth .le. -1. ) omega = 0.
    if ( tanLatMonth .ge. 1. )  omega = pi
    N_up = 24.*omega / pi
    K_up(month) = N_up/12. * days(month-1)/30.
  enddo
  !! Compute potential evapotranspiration series
  month = 1
  do i = 0, n - 1
    if ( (tempSeries(i) .gt. 0.) .and. (J .gt. 0.) ) then
      etpSeries(i) = K_up(month) * 16. * (10*tempSeries(i)/J)**c
    else
      etpSeries(i) = 0.
    endif

    if (month .lt. 12) then
      month = month + 1
    else
      month = 1
    endif
  enddo

End Subroutine thornthwaite

!-------------------------------------------------------------------
Subroutine Hargreaves (Tx, Tm, n, lat, etpSeries)
  ! ETo = 0.0023 * RA * (Tave + 17.8) * TD**0.5
  ! ETo = Krs * RA * (Tave + HT) * TD**HE
  Integer, Parameter :: Nmon=12, Nlat=51
  Real, Parameter :: Krs=0.0023, HT=17.8, HE=0.5
  real, parameter :: days (Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/)
  integer :: n
  real :: Tx(0:n-1), Tm(0:n-1)
  real :: Tave (0:n-1), TD(0:n-1), lat, etpSeries (0:n-1)
  Real :: RA0(Nlat,Nmon), Lat0(Nlat), RA(Nmon)
  integer :: month, i, k

  DATA Lat0 &
  /-50.0,-48.0,-46.0,-44.0,-42.0,-40.0,-38.0,-36.0,-34.0,-32.0,-30.0,-28.0,-26.0,-24.0,-22.0,-20.0,-18.0,-16.0,-14.0,-12.0,-10.0, -8.0, -6.0, -4.0, -2.0,  0.0, &
     2.0,  4.0,  6.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 22.0, 24.0, 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 38.0, 40.0, 42.0, 44.0, 46.0, 48.0, 50.0/
  DATA RA0 &
  /17.5, 17.6, 17.7, 17.8, 17.8, 17.9, 17.9, 17.9, 17.8, 17.8, 17.8, 17.7, 17.6, 17.5, 17.4, 17.3, 17.1, 16.9, 16.7, 16.6, 16.4, 16.1, 15.8, 15.5, 15.3, 15.0, &
   14.7, 14.3, 13.9, 13.6, 13.2, 12.8, 12.4, 12.0, 11.6, 11.2, 10.7, 10.2,  9.8,  9.3,  8.8,  8.3,  7.9,  7.4,  6.9,  6.4,  5.9,  5.3,  4.9,  4.3,  3.8, 14.7, &
   14.9, 15.1, 15.3, 15.5, 15.7, 15.8, 16.0, 16.1, 16.2, 16.4, 16.4, 16.4, 16.5, 16.5, 16.5, 16.5, 16.4, 16.4, 16.3, 16.3, 16.1, 16.0, 15.8, 15.7, 15.5, 15.3, &
   15.0, 14.8, 14.5, 14.2, 13.9, 13.6, 13.3, 13.0, 12.7, 12.3, 11.9, 11.5, 11.1, 10.7, 10.2,  9.8,  9.4,  9.0,  8.6,  8.1,  7.6,  7.1,  6.6,  6.1, 10.9, 11.2, &
   11.5, 11.9, 12.2, 12.5, 12.8, 13.2, 13.5, 13.8, 14.0, 14.3, 14.4, 14.6, 14.8, 15.0, 15.1, 15.2, 15.3, 15.4, 15.5, 15.5, 15.6, 15.6, 15.7, 15.7, 15.6, 15.5, &
   15.4, 15.3, 15.3, 15.1, 14.9, 14.7, 14.6, 14.4, 14.2, 13.9, 13.7, 13.4, 13.1, 12.8, 12.4, 12.1, 11.8, 11.4, 11.0, 10.6, 10.2,  9.8,  9.4,  7.0,  7.5,  7.9, &
    8.4,  8.8,  9.2,  9.6, 10.1, 10.5, 10.9, 11.3, 11.6, 12.0, 12.3, 12.6, 13.0, 13.2, 13.5, 13.7, 14.0, 14.2, 14.4, 14.7, 14.9, 15.1, 15.3, 15.3, 15.5, 15.4, &
   15.6, 15.7, 15.7, 15.7, 15.6, 15.6, 15.6, 15.5, 15.4, 15.3, 15.3, 15.2, 15.0, 14.8, 14.7, 14.5, 14.3, 14.0, 13.7, 13.3, 13.0, 12.7,  4.2,  4.7,  5.2,  5.7, &
    6.1,  6.6,  7.1,  7.5,  8.0,  8.5,  8.9,  9.3,  9.7, 10.2, 10.6, 11.0, 11.4, 11.7, 12.1, 12.5, 12.8, 13.1, 13.4, 13.8, 14.1, 14.4, 14.6, 14.9, 15.1, 15.3, &
   15.5, 15.7, 15.8, 16.0, 16.1, 16.3, 16.3, 16.4, 16.4, 16.5, 16.5, 16.5, 16.5, 16.4, 16.4, 16.4, 16.2, 16.1, 16.0, 15.9, 15.8,  3.1,  3.5,  4.0,  4.4,  4.9, &
    5.3,  5.8,  6.3,  6.8,  7.3,  7.8,  8.2,  8.7,  9.1,  9.6, 10.0, 10.4, 10.8, 11.2, 11.6, 12.0, 12.4, 12.8, 13.2, 13.5, 13.9, 14.2, 14.4, 14.7, 15.0, 15.3, &
   15.5, 15.7, 15.9, 16.1, 16.4, 16.4, 16.6, 16.7, 16.8, 17.0, 17.0, 17.1, 17.2, 17.2, 17.3, 17.3, 17.2, 17.2, 17.2, 17.1,  3.5,  4.0,  4.4,  4.9,  5.4,  5.9, &
    6.3,  6.8,  7.2,  7.7,  8.1,  8.6,  9.1,  9.5, 10.0, 10.4, 10.8, 11.2, 11.6, 12.0, 12.4, 12.7, 13.1, 13.4, 13.7, 14.1, 14.3, 14.6, 14.9, 15.1, 15.3, 15.5, &
   15.7, 15.9, 16.1, 16.3, 16.4, 16.5, 16.6, 16.7, 16.8, 16.8, 16.8, 16.7, 16.7, 16.7, 16.7, 16.6, 16.6, 16.5, 16.4,  5.5,  6.0,  6.5,  6.9,  7.4,  7.9,  8.3, &
    8.8,  9.2,  9.6, 10.1, 10.4, 10.9, 11.2, 11.6, 12.0, 12.3, 12.6, 12.9, 13.2, 13.5, 13.7, 14.0, 14.3, 14.5, 14.8, 14.9, 15.1, 15.2, 15.4, 15.5, 15.6, 15.7, &
   15.7, 15.8, 15.9, 15.8, 15.8, 15.7, 15.7, 15.7, 15.6, 15.5, 15.4, 15.3, 15.2, 15.0, 14.7, 14.5, 14.3, 14.1,  8.9,  9.3,  9.7, 10.2, 10.6, 11.0, 11.4, 11.7, &
   12.0, 12.4, 12.7, 13.0, 13.2, 13.4, 13.7, 13.9, 14.1, 14.3, 14.5, 14.7, 14.8, 14.9, 15.0, 15.1, 15.2, 15.3, 15.3, 15.3, 15.3, 15.3, 15.3, 15.2, 15.1, 15.0, &
   14.9, 14.8, 14.6, 14.5, 14.3, 14.1, 13.9, 13.6, 13.4, 13.1, 12.8, 12.5, 12.2, 11.9, 11.5, 11.2, 10.9, 12.9, 13.2, 13.4, 13.7, 14.0, 14.2, 14.4, 14.6, 14.9, &
   15.1, 15.3, 15.4, 15.5, 15.6, 15.7, 15.8, 15.8, 15.8, 15.8, 15.8, 15.9, 15.8, 15.7, 15.6, 15.5, 15.4, 15.3, 15.1, 15.0, 14.8, 14.7, 14.4, 14.1, 13.9, 13.6, &
   13.3, 13.0, 12.6, 12.3, 12.0, 11.6, 11.2, 10.8, 10.6, 10.0,  9.6,  9.1,  8.7,  8.3,  7.8,  7.4, 16.5, 16.6, 16.7, 16.7, 16.8, 16.9, 17.0, 17.0, 17.1, 17.2, &
   17.3, 17.2, 17.2, 17.1, 17.0, 17.0, 16.8, 16.7, 16.5, 16.4, 16.2, 16.0, 15.8, 15.5, 15.3, 15.1, 14.8, 14.5, 14.2, 13.9, 13.6, 13.3, 12.8, 12.4, 12.0, 11.6, &
   11.1, 10.7, 10.3,  9.9,  9.5,  9.0,  8.5,  8.0,  7.5,  7.0,  6.5,  6.0,  5.5,  5.0,  4.5, 18.2, 18.2, 18.3, 18.3, 18.3, 18.3, 18.3, 18.2, 18.2, 18.1, 18.1, &
   17.9, 17.8, 17.7, 17.5, 17.4, 17.1, 16.8, 16.6, 16.5, 16.2, 16.0, 15.7, 15.4, 15.1, 14.8, 14.4, 14.1, 13.7, 13.3, 12.9, 12.5, 12.0, 11.6, 11.1, 10.7, 10.2, &
    9.7,  9.3,  8.8,  8.3,  7.8,  7.2,  6.6,  6.1,  5.7,  5.2,  4.7,  4.3,  3.7,  3.2 /
   !-------------------------------------------------
   ! Calculate RA = f(lat,mon)
   Do k=1,Nmon
      Do i=1,Nlat-1
         if (lat>Lat0(i).AND.lat<=Lat0(i+1)) exit
      Enddo
      RA(k)=RA0(i,k)+(RA0(i+1,k)-RA0(i,k))/(Lat0(i+1)-Lat0(i))*(lat-Lat0(i))
   Enddo
   ! Calculate PET
   k=0
   Do i=0,n-1
      if (Tx(i)>0.AND.Tm(i)>0) then
         Tave(i)=(Tx(i)+Tm(i))/2.0; TD(i)=Tx(i)-Tm(i)
      endif
      k=mod(i+1,Nmon); if (k==0) k=Nmon
      etpSeries(i) = Krs * RA(k) * (Tave(i) + HT) * TD(i)**HE  ! mm/day
      etpSeries(i) = etpSeries(i) * days(k)                    ! mm/month
   Enddo
End Subroutine
!-------------------------------------------------------------------