!------------------------------------------------------------------------------ Subroutine Cal_SPEI Integer :: ist, i, acumRegistros, TimeScale Real, Dimension (0:Ntime-1) :: balanceSeries, acumSeries, speiSeries Loop_1: Do ist=1, Nst Rain(0:Ntime-1)=DatBL(1:Ntime,ist); T2m(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(1:Ntime) ! pause Do i=1,NLen 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) spi(ist,TimeScale:Ntime,i) = speiSeries(0:acumRegistros-1) 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 Cal_SPEI !---------------------------------------------------------------------------------! 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 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 !------------------------------------------------------------------------------ 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 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 !---------------------------------------------------------------------------------! !! 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: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 !---------------------------------------------------------------------------------! 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 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 !---------------------------------------------------------------------------------! !! 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