! Tinh SPEI cho mot monthly timeseries of Rainfall and Temp cua mot station ! So lieu dau vao bat buoc la: ! + Ntime, Nmon, TimeScale, T2m(Ntime), Rain(Ntime), Lat, use_baseline, Miss ! + Neu use_baseline=.true. thi can SL cua logLogisticParams(0:Nmon,0:2) ! =.false. thi logLogisticParams(0:Nmon,0:2) se la dau ra ! Outputs: ! + SPEI_out(Ntime) ! + logLogisticParams(0:Nmon,0:2) neu use_baseline=.false. !---------------------------------------------------------------------------------! !Subroutine Calc_SPEI_Thorn_Daily_One_Stn (T2m, Rain, Lat, Ntime, Nmon, TimeScale, logLogisticParams, use_baseline, SPEI_out, Miss) Subroutine Calc_SPEI_Thorn_Daily_One_Stn (T2m, Rain, Lat, Ntime, TimeScale, SPEI_out, Miss) ! tempSeries: Timeseries of Temperature (T2m) [Ntime] ! rainSeries: Timeseries of Rainfall (Rain) [Ntime] ! lat: Latitude of station ! numRegistros: Timeseries length (Ntime) ! acumRegistros: Number of SPEI's valid values (=Ntime-TimeScale = Nspei) ! acumulated: TimeScale ! speiSeries: SPEI values [Ntime] ! ==> Call Calc_SPEI_One_Stn (T2m, Rain, Lat, Ntime, TimeScale, speiSeries) <====== Integer :: Ntime, TimeScale, Nmon, Nt1, Nspei integer :: acumRegistros, indice, jndice real :: lat, Miss, sumx Real, Dimension (1:Ntime) :: T2m, Rain, SPEI_out !Real, Dimension (0:Ntime-1) :: rainSeries, tempSeries, etpSeries, balanceSeries, acumSeries, speiSeries Real, Dimension (0:Ntime-1) :: etpSeries, balanceSeries, acumSeries, speiSeries !Real :: logLogisticParams(0:Nmon,0:2) !Logical :: use_baseline, L1, L2 !---------------------------------------------------------------------------------------- ! L1 = Empty_Arr1_spei(Rain, Ntime, Miss); L2 = Empty_Arr1_spei(T2m, Ntime, Miss) !numRegistros = Ntime; !acumulated = TimeScale ! Su dung ky hieu cua original version SPEI_out(:) = Miss; !if (Empty_Arr1_spei(Rain, Ntime, Miss).OR.Empty_Arr1_spei(T2m, Ntime, Miss)) then !if (L1.OR.L2) then ! SPEI_out(:) = Miss; return !endif !tempSeries(0:Ntime-1) = T2m (1:Ntime) !rainSeries(0:Ntime-1) = Rain(1:Ntime) !! Compute the climatic balance: precipitation minus potential evapotranspiration !if ( (tempSeries(1) .ne. 0.) .and. (tempSeries(2) .ne. 0.) ) then !print*,"In SPEI Sub: TimeScale ",tempSeries(1),tempSeries(2) !if ( (tempSeries(1) .ne. Miss) .and. (tempSeries(2) .ne. Miss) ) then !if (.not.L2) then ! Call thornthwaite (tempSeries, Ntime, lat, etpSeries, Miss) Call Thornthwaite_Daily (T2m, Ntime, lat, etpSeries, Miss) do indice = 0, Ntime - 1 if (Rain(indice+1)/=Miss.AND.etpSeries(indice)/=Miss) then ! TP-2023 balanceSeries(indice) = Rain(indice+1) - etpSeries(indice) else balanceSeries(indice) = Miss endif !print*,Rain(indice+1),etpSeries(indice),balanceSeries(indice) enddo !pause !else ! do indice = 0, Ntime - 1 ! balanceSeries(indice) = rainSeries(indice) ! enddo !endif !---------------------------------------------------------------------------------------- !! Compute the cumulative series !print*,"In SPEI Sub: TimeScale ",TimeScale,acumulated acumSeries(:) = Miss ! 0.; Nt1 = TimeScale do indice = TimeScale - 1, Ntime -1 ! Accumulate water balance over TimeScale sumx = 0.0 do jndice = 0, TimeScale - 1 if (balanceSeries (indice - jndice)/=Miss) then ! TP-2023 !acumSeries(indice - TimeScale + 1) = & ! Indix value from 0 !acumSeries(indice - TimeScale + 1) + balanceSeries (indice - jndice) sumx = sumx + balanceSeries (indice - jndice) endif enddo acumSeries(indice) = sumx ! acumSeries(0:TimeScale-2) = Miss enddo !---------------------------------------------------------------------------------------- !! Compute the SPEI series Call spei(acumSeries, Ntime, speiSeries, Miss) !, logLogisticParams, use_baseline) SPEI_out(1:Ntime) = speiSeries(0:Ntime-1) !SPEI_out(Nt1:Ntime) = speiSeries(0:Nspei-1) ! print*,"Ntime,SPEI: ",Ntime,acumRegistros,Nt1,SPEI_out(Ntime) ! print"(24F6.3)",SPEI_out(1:Ntime) ! print"(12F10.3)", ((logLogisticParams(j,k),j=1,Nmon),k=0,2) !pause End Subroutine Calc_SPEI_Thorn_Daily_One_Stn !------------------------------------------------------------------------------ !Subroutine Calc_SPEI_Harg_One_Stn (Tmx, Tmn, Rain, Lat, Ntime, Nmon, TimeScale, logLogisticParams, use_baseline, SPEI_out, Miss) Subroutine Calc_SPEI_Harg_Daily_One_Stn (Tmx, Tmn, Rain, Lat, Ntime, TimeScale, SPEI_out, Miss) Integer :: Ntime, TimeScale integer :: acumulated, numRegistros, acumRegistros, indice, jndice, Ntime1 Integer :: Nspei, Nt1 real :: lat, Miss, sumx real, dimension (1:Ntime) :: Tmx, Tmn, Rain, SPEI_out !real, dimension (0:Ntime-1) :: rainSeries, Tx, Tm, etpSeries, balanceSeries, acumSeries, speiSeries real, dimension (0:Ntime-1) :: etpSeries, balanceSeries, acumSeries, speiSeries !Logical :: use_baseline, L1, L2, L3 !SPEI_out(:) = Miss; !L1 = Empty_Arr1_spei(Rain, Ntime, Miss) !L2 = Empty_Arr1_spei(Tmx, Ntime, Miss) !L3 = Empty_Arr1_spei(Tmn, Ntime, Miss) !if (Empty_Arr1_spei(Rain, Ntime, Miss).OR.Empty_Arr1_spei(Tmx, Ntime, Miss).OR.Empty_Arr1_spei(Tmn, Ntime, Miss)) then !if (L1.OR.L2.OR.L3) then ! SPEI_out(:) = Miss; return !endif !Tx(0:Ntime-1) = Tmx(1:Ntime) !Tm(0:Ntime-1) = Tmn(1:Ntime) !rainSeries(0:Ntime-1) = Rain(1:Ntime) !Ntime1 = Ntime -1 !numRegistros = Ntime; acumulated = TimeScale SPEI_out(:) = Miss; !! Compute the climatic balance: precipitation minus potential evapotranspiration !if (Tx(1)/=Miss.AND.Tx(2)/=0.AND.Tm(1)/=0.AND.Tm(2)/=0) then !if (.not.Empty_Arr1_spei (Tmx, Ntime, Miss).AND..not.Empty_Arr1_spei (Tmn, Ntime, Miss)) then !if (.not.L2.AND..not.L3) then ! Call Hargreaves (Tx, Tm, numRegistros, lat, etpSeries) Call Hargreaves_daily (Tmx, Tmn, Ntime, lat, etpSeries, Miss) do indice = 0, Ntime - 1 if (Rain(indice+1)/=Miss.AND.etpSeries(indice)/=Miss) then balanceSeries(indice) = Rain(indice+1) - etpSeries(indice) else balanceSeries(indice) = Miss endif enddo !else ! do indice = 0, numRegistros - 1 ! balanceSeries(indice) = rainSeries(indice) ! enddo !endif !---------------------------------------------------------------------------------------- !! Compute the cumulative series ! acumRegistros = numRegistros - acumulated + 1; acumSeries(:) = 0. ! Nspei = numRegistros - acumulated + 1; acumSeries(:) = 0. acumSeries(:) = Miss !; Nt1 = TimeScale do indice = TimeScale - 1, Ntime -1 sumx = 0 do jndice = 0, TimeScale - 1 if (balanceSeries (indice - jndice)/=Miss) then ! TP-2023 !acumSeries(indice - TimeScale + 1) = & ! Indix value from 0 !acumSeries(indice - TimeScale + 1) + balanceSeries (indice - jndice) sumx = sumx + balanceSeries (indice - jndice) endif enddo acumSeries(indice) = sumx ! acumSeries(0:TimeScale-2) = Miss enddo !---------------------------------------------------------------------------------------- !! Compute the SPEI series !Call spei(acumSeries, acumRegistros, NMonth, speiSeries) !print"(20F7.1)",acumSeries(0:Ntime-1);print*,"Acumulated BW"; pause Call spei(acumSeries, Ntime, speiSeries, Miss) !, logLogisticParams, use_baseline) SPEI_out(1:Ntime) = speiSeries(0:Ntime-1) !SPEI_out(Nt1:Ntime) = speiSeries(0:Nspei-1) End Subroutine Calc_SPEI_Harg_Daily_One_Stn !------------------------------------------------------------------------------ Logical Function Empty_Arr1_spei (X, N, Miss) ! Tra ve .True. neu X(N) khuyet toan bo SL Integer :: N, i, k Real :: X(1:N), Miss k = 0 ! Dem so phan tu khuyet cua X Do i=1, N if (X(i) == Miss) k = k + 1 Enddo if (k > N/5) then ; Empty_Arr1_spei = .True. ; else ; Empty_Arr1_spei = .False. ; endif End Function Empty_Arr1_spei !------------------------------------------------------------------------------------------ Subroutine Monthly_Data_Filter (X,Ntime,Nyear,Y,Ny,imon,Mean,Miss) Integer, Parameter :: N365=365, Nmon=12 Integer, parameter :: DayOfMn(Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/) Integer :: AccDays !(Nmon) Integer :: Ntime, Nyear, Ny, it, iday, imon, i, iyr Real :: X(0:Ntime-1), Y(0:Ntime-1), Miss Real :: tmp(Nmon), ic, HI, day Real, Allocatable :: Mean(:) Nyear = Ntime/N365; if (Allocated (Mean)) Deallocate (Mean) Allocate (Mean(Nyear)) ! Monthly values AccDays = 0 Do i=1,imon AccDays = AccDays + DayOfMn(i) ! So thu tu ngay cuoi thang trong nam (Jan=31, Feb=59, etc) Enddo Ny = 0; iyr = 1; Mean(:) = 0 Do it=1,Ntime iday=mod(it,N365); if (iday==0) then; iday=N365; iyr = iyr+1; endif ! Tang 1 nam if (imon==1) then if (iday<=AccDays) then if (X(it-1)/=Miss) then Ny=Ny+1; Y(Ny-1)=X(it-1) Mean(iyr) = Mean(iyr)+X(it-1) ! Sum values in month endif endif else if (iday>AccDays-DayOfMn(imon).and.iday<=AccDays) then if (X(it-1)/=Miss) then Ny=Ny+1; Y(Ny-1)=X(it-1) Mean(iyr) = Mean(iyr)+X(it-1) ! Sum values in month endif endif endif Enddo End Subroutine Monthly_Data_Filter !---------------------------------------------------------------------------------! Subroutine spei(dataSeries, N, speiSeries, Miss) !, logLogisticParams, use_baseline) integer :: N real, dimension (0:n-1) :: dataSeries, speiSeries Integer, Parameter :: N365=365, Nmon=12 Integer, parameter :: DayOfMn(Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/) integer :: i, j, k, it, nSeason Integer :: AccDays, iday, Nyear Real, Allocatable :: Mean(:) real :: seasonSeries (0:n-1), beta (0:2), logLogisticParams(0:Nmon,0:2) Real :: Miss !------------------------------------------------------------------------------ speiSeries(:) = Miss do j = 1, Nmon !seasons !k = 0 !do i = j - 1, n - 1, seasons ! seasonSeries (k) = dataSeries (i) ! k = k + 1 !enddo !nSeason = k ! nSeason = k - 1 ! TP them -1 Call Monthly_Data_Filter (dataSeries,N,Nyear,seasonSeries,nSeason,j,Mean,Miss) ! nSeason = So gia tri (daily) cua seasonSeries !!! Call upward (seasonSeries, nSeason) ! Sap xep chuoi tang dan (Daily data) Call upward (Mean,Nyear) ! Sap xep chuoi tang dan (Monthly data) !! Compute probability weighted moments !!! Call pwm (seasonSeries, nSeason, beta, 0., 0., 0) ! Daily data Call pwm (Mean,Nyear,beta, 0., 0., 0) ! Monthly data !!print*,"beta",beta(0), beta(1), beta(2) !! Fit a Log Logistic probability function Call logLogisticFit (beta, logLogisticParams(j,:)) !print*,"logLogisticFit:",j, logLogisticParams(j,0), logLogisticParams(j,1), logLogisticParams(j,2) !! Calculate the standardized values ! do i = j - 1, n - 1, Nmon ! seasons ! speiSeries (i) = logLogisticCDF (dataSeries (i), logLogisticParams (j,:)) ! !print*,standardGaussianInvCDF (speiSeries(i)) ! speiSeries (i) = -standardGaussianInvCDF (speiSeries(i)) ! !print*,"CHECK: ",j,n-1,i,speiSeries (i) ! enddo AccDays = 0 Do i=1, j ! Month j AccDays = AccDays + DayOfMn(i) ! So thu tu ngay cuoi thang trong nam (Jan=31, Feb=59, etc) Enddo Do it=1,N ! Apply distribution of monthly data to daily data iday=mod(it,N365); if (iday==0) iday=N365 if (j==1) then ! month j=1 if (iday<=AccDays) then if (dataSeries(it-1)/=Miss) then speiSeries(it-1) = logLogisticCDF (dataSeries (it-1), logLogisticParams (j,:)) speiSeries(it-1) = -standardGaussianInvCDF (speiSeries(it-1)) endif endif else if (iday>AccDays-DayOfMn(j).and.iday<=AccDays) then if (dataSeries(it-1)/=Miss) then speiSeries(it-1) = logLogisticCDF (dataSeries (it-1), logLogisticParams (j,:)) speiSeries(it-1) = -standardGaussianInvCDF (speiSeries(it-1)) endif endif endif 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) ! 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=Tm(it)) then etpSeries(it-1) = PET_daily_Harg (Tx(it), Tm(it), lat, iday) else if (Tx(it)AccDays(imon-1).and.iday<=AccDays(imon)) then if (T2m(it)/=Miss) then ic=ic+1; tmp(imon)=tmp(imon)+T2m(it) endif endif endif Enddo tmp(imon) = tmp(imon)/ic ! Climatological Monthly Mean !print*,tmp(imon) Enddo !pause HI = 0.0 Do imon=1,Nmon ic = Max (0.0, tmp(imon)) HI = HI + (ic/5.0) ** 1.514 Enddo Heat_Index = HI End Function Heat_Index !------------------------------------------------------------------------------------------ Subroutine Thornthwaite_Daily (T2m, Ntime, lat, etpSeries, Miss) ! Calculate Daily PET using Thornthwaite method ! Input: T2m(Ntime) daily mean temp in C degree, lat (in degree) ! Output: etpSeries(0:Ntime-1) daily potential evapolation ! Note: Ntime must begin counting from 1st January Integer, Parameter :: N365=365, Nmon=12 Integer, parameter :: DayOfMn(Nmon) = (/31,28,31,30,31,30,31,31,30,31,30,31/) Integer :: AccDays(Nmon) Integer :: Ntime, it, iday, imon Real :: T2m(Ntime), etpSeries(0:Ntime-1), lat, Miss Real :: tmp(Nmon), ic, HI, day HI = Heat_Index (T2m, Ntime, Miss) Do it=1,Ntime iday = mod(it,N365); if (iday==0) iday = N365 if (T2m(it)/=Miss) then etpSeries(it-1) = PET_daily_Thorn (T2m(it), HI, lat, iday) else etpSeries(it-1) = Miss endif Enddo End Subroutine Thornthwaite_Daily !------------------------------------------------------------------------------------------ Real Function PET_daily_Thorn (t2m, I, lat, iday) ! mm/day ! Calculate Potential Evaporation for a certain day of the year (day: 1-365) at latitude (lat in degree) ! t2m is daily mean temperature ! I is heat index calculated from monthly climatological mean temperature: I = SUM from Jan to Dec of (max[0, Tmean]/5**1.514) ! Input: T2m (of iday), lat (degree), iday (day of the year), I (Heat Index) Integer :: iday Real :: t2m, I, lat, day, pet, tanlatday Real :: a, N, ws, phi, delta, pi, pi2d365, I2, I3, C, C1, C2, T2 Real, parameter :: a1=0.000000675, a2=0.0000771, a3=0.01792, a4=0.49239 pi = 3.141592654 phi = pi/180.0*lat ! Latitude in radian pi2d365 = 2.0*pi/365. day = iday delta = 0.409*sin(pi2d365*day-1.39) ! solar declination angle in radian (xich vi) tanlatday = tan(phi)*tan(delta) if ( (tanlatday > -1.) .and. (tanlatday < 1.) ) ws = acos(-tanlatday) if ( tanlatday <= -1. ) ws = 0. if ( tanlatday >= 1. ) ws = pi ! if (abs(delta)-pi/2. < 1.0E-6) then ! ws = 0.1 ! else ! ws = acos(-tan(phi)*tan(delta)) ! sunset hour angle, sunset ! endif N = 24./pi*ws ! daylight hours C = N/360.0 ! Convert to daily I2 = I+I; I3=I*I*I !a = a1*I3 - a2*I2 + a3+I + a4 a = 0.000000675*I3 - 0.0000771*I2 + 0.01792*I + 0.49239 C1 = (10.0*t2m/I)**a ! In case t2m in between 0 - 26 T2 = t2m*t2m C2 = -415.85 + 32.24*t2m - 0.43*T2 ! In case t2m > 26 !print"(A,10F12.3)","==== a,C,c1,c2,I ====>",a,C,c1,c2,I if (t2m < 0) then pet = 0. else if (t2m <= 26.0) then pet = 16.0 * C * C1 else pet = C * C2 endif PET_daily_Thorn = pet End Function PET_daily_Thorn !------------------------------------------------------------------------------------------ Real Function PET_daily_Harg (tx, tm, lat, iday) ! mm/day ! Calculate Potential Evaporation for a certain day of the year (day: 1-365) at latitude (lat in degree) ! tx, tm are daily maximum and minimum temperature, respectively ! https://github.com/jernejjevsenak/calculate-daily-SPEI/blob/main/HargreavesPET_daily.R ! Input: Tx, Tm (of iday), lat (degree), iday (integer) day of the year Integer :: iday ! day of the year (1-365) Real :: tx, tm, lat, day Real :: tavg, tdif, pi, phi, delta, Ra, dr, ws, pet, pi2d365, tanlatday pi = 3.141592654 phi = pi/180.0*lat ! Latitude in radian tdif = tx - tm ! range temp tavg = (tx + tm)/2.0 ! mean temp pi2d365 = 2.0*pi/365. ! day = iday delta = 0.409*sin(pi2d365*day-1.39) ! solar declination angle in radian (xich vi) ! https://www.pveducation.org/pvcdrom/properties-of-sunlight/declination-angle dr = 1.0 + 0.033*cos(pi2d365*day) ! relative distance Earth-Sun (1 astronomical unit - AU) tanlatday = tan(phi)*tan(delta) if ( (tanlatday > -1.) .and. (tanlatday < 1.) ) ws = acos(-tanlatday) if ( tanlatday <= -1. ) ws = 0. if ( tanlatday >= 1. ) ws = pi !if (abs(delta)-pi/2. < 1.0E-6) then ! ws = 0.1 !else ! ws = acos(-tan(phi)*tan(delta)) ! sunset hour angle, sunset !endif Ra = (24*60/pi)*0.0820*dr*((ws*sin(phi)*sin(delta))+(cos(phi)*cos(delta)*sin(ws))) if (Ra < 0) Ra = 0.0 pet = 0.0023*(tavg + 17.8)*sqrt(tdif) * Ra !if (IsNaN(pet)) print*,ws,ra,tdif if (pet < 0) then PET_daily_Harg = 0.0 else PET_daily_Harg = pet endif End Function PET_daily_Harg !------------------------------------------------------------------------------