thornthwaite
Estimate the potential evapotranspiration (PET) via the Thornthwaite method.
Available in version 6.3.0 and later.
Prototype
function thornthwaite ( t : numeric, lat : numeric, opt [1] : logical, dim [1] : integer ) return_val : float or double
Arguments
tMonthly near surface temperature (C). The number of months (ntim) must be a multiple of 12. Further, it is assumed that the data start in January. Common dimensionalities include: t(ntim), t(ntim,ncol), t(ntim,nlat,mlon), t(ncase,ntim), t(ncase,ntim,nlat,mlon), t(ncase,ntim,klon).
latLatitude(s) of the current station or grid point(s).
- t(ntim), then lat is a scalar and dim=0
- t(ntim,ncol), and the grid is unstructured (eg: spectral element); then lat(ncol) and dim=0
- t(ntim,nlat,mlon), and the grid is rectilinear, then lat(nlat) and dim=0
- t(ntim,nlat,mlon), and the grid is curvilinear, then lat(nlat,mlon) and dim=0
- t(ncase,ntim), then lat is a scalar and dim=1
- t(ncase,ntim,ncol), and the grid is unstructured (eg: spectral element); then lat(ncol) and dim=1
- t(ncase,ntim,nlat,mlon), and the grid is curvilinear, then lat(nlat,mlon) and dim=1
- t(ncase,ntim,nlat,mlon), and the grid is curvilinear, then lat(nlat,mlon) and dim=1
Not used. Set to False.
dimThe dimension of t to be used to estimate the potential evapotranspiration (PET). This is the time dimension.
Return value
The returned PET will be the same shape, size and type as t. The units are millimeters/month (ie: monthly total).
Description
Potential evapotranspiration (PET) is estimated via the method of Thornthwaite. The only variable used is monthly temperature [t (C)]. The calculated climatological means are used to calculate an annual temperature efficiency index. The latitude is used to adjust the number of sunlight hours over the course of the year.
This method has been the subject of many criticisms. Generally, the Thornthwaite method underestimates PET in the arid areas and over estimates PET in the humid areas.
Also, there is confusion about the difference between potential evapotranspiration (PET) and reference evapotranspiration (ETo). Irmak and Haman (2014) provide a nice description in their paper entitled Evapotranspiration: Potential or Reference?. [UF: IFAS Extension: ABE 343].
References:
C.W. Thornthwaite (1948):
An approach toward a rational classification of climate, Geographical Review 38 (1): 55-94.
http://dx.doi.org/10.2307/21073
A.R. Pereira, Â. Paes De Camargo (1989):
An analysis of the criticism of thornthwaite's equation for estimating potential evapotranspiration
Agricultural and Forest Meteorology, Volume 46, Issue 1, Pages 149-157
http://dx.doi.org/10.1016/0168-1923(89)90118-4
F. Hashemi, M.T. Habibian (1979):
Limitations of temperature-based methods in estimating crop evapotranspiration
in arid-zone agricultural development projects
Agricultural Meteorology, Volume 20, Issue 3, June 1979, Pages 237-247
http://dx.doi.org/10.1016/0002-1571(79)90025-6
See Also
dim_spi_n, dim_spi3_n, dim_gamfit_n
Examples
Example 1
Sellers (1969), Physical Climatology, pg 173 includes an example for one year at Aspendale, Australia. Compute potential evapotranspiration (mm/month) and compare with Sellers Table 23 which are in units of mm/day.
tmp = (/23.3, 21.1, 19.6, 17.2, 12.6, 10.9 \
,10.0, 11.0, 13.0, 15.8, 17.8, 20.1/)
tmp@units = "degC"
tmp@long_name = "Sellers, Physical Climatology; 1969, p173"
printVarSummary(tmp)
print("min="+min(tmp)+" max="+max(tmp))
latAsp = -38.0 ; Latitude of Aspendale, Australia
pet = thornthwaite(tmp, latAsp, False, 0)
pet@long_name = "PET: Aspendale, Australia"
pet@units = "mm/month" ; ie, monthly total
; Compare with Sellers Table 23, pg 173 which are in units of mm/day
dymo = (/31,28,31,30,31,30,31,31,30,31,30,31/)
pet = pet/dymo ; mm/month => mm/day
pet@units = "mm/day"
pet_sellers = (/4.47, 3.51, 2.80, 2.04, 1.09, 0.80 \
,0.74, 0.89, 1.50, 2.09, 2.73, 3.47 /)
pet_sellers@long_name = "PET: Sellers Table 23"
pet_sellers@units = "mm/day"
printVarSummary(pet_sellers)
print("min="+min(pet_sellers)+" max="+max(pet_sellers))
print(" ")
print(pet_sellers+" "+sprintf("%6.2f",pet) )
The (edited) output looks like:
Variable: pet_sellers
Type: float
Total Size: 48 bytes
12 values
Number of Dimensions: 1
Dimensions and sizes: [12]
Coordinates:
Number Of Attributes: 2
units : mm/day
long_name : PET: Sellers Table 23
(0) min=0.74 max=4.47
Sellers NCL
(0) 4.47 4.32
(1) 3.51 3.46
(2) 2.80 2.81
(3) 2.04 2.03
(4) 1.09 1.11
(5) 0.80 0.83
(6) 0.74 0.74
(7) 0.89 0.94
(8) 1.50 1.37
(9) 2.09 2.09
(10) 2.73 2.73
(11) 3.47 3.45
Example 2
Read an ascii file containing Palmer's Table 1 ("Hydrologic accounting for central Iowa"). There is no indication what latitude Palmer used for 'central Iowa.' This example uses 42N. This is from Palmer's (1965) classic note: Meteorological Drought: https://www.ncdc.noaa.gov/temp-and-precip/drought/docs/palmer.pdf
diri = "./" fili = "palmer_table1.txt" ; Palmer Table 1: yyyymm, T(F), P(in), PE(in) data = readAsciiTable(diri+fili, 4, "float", 3) yyyymm = toint(data(:,0)) tf = data(:,1) ; Farenheit prcin = data(:,2) ; Inches petin = data(:,3) ; Inches lat = 42.0 ; central Iowa (approximate) ; convert 'tf' to degC [required by 'thornthwaite'] PETncl = thornthwaite((tf-32)*0.555556, lat, False 0) ; mm/month PETncl = PETncl*0.03937008 ; mm => inches for comparison with Palmer values petdif = PETncl-petin ; difference for print print(yyyymm+ sprintf("%8.2f", tf)+\ sprintf("%8.2f", prcin)+\ sprintf("%8.2f", petin)+\ sprintf("%8.2f",PETncl)+\ sprintf("%8.2f",petdif)) The (slightly edited) output is as follows. All precipitation rekated values are in inches for comparison.yyyymm Tf Pin PE PETncl difference (0) 195301 32.30 0.90 0.01 0.01 -0.00 (1) 195302 22.60 0.21 0.00 0.00 0.00 (2) 195303 36.30 3.22 0.31 0.23 -0.08 (3) 195304 48.50 1.16 1.63 1.45 -0.18 (4) 195305 60.60 5.60 3.56 3.36 -0.20 (5) 195306 77.80 1.03 6.37 6.25 -0.12 (6) 195307 76.30 3.57 6.22 6.09 -0.13 (7) 195308 70.30 1.84 4.84 4.70 -0.14 (8) 195309 69.70 3.57 4.18 4.01 -0.17 (9) 195310 49.80 1.95 1.52 1.37 -0.15 (10) 195311 37.90 0.26 0.34 0.28 -0.06 (11) 195312 26.90 0.89 0.00 0.00 0.00 (12) 195401 26.10 1.34 0.00 0.00 0.00 (13) 195402 25.10 0.59 0.00 0.00 0.00 (14) 195403 34.40 1.01 0.09 0.11 0.02 (15) 195404 50.10 0.61 1.82 1.64 -0.18 (16) 195405 70.00 0.76 5.08 4.87 -0.21 (17) 195406 78.40 2.10 6.53 6.36 -0.17 (18) 195407 79.30 4.68 6.82 6.63 -0.19 (19) 195408 72.90 2.83 5.25 5.12 -0.13 (20) 195409 61.10 5.59 3.02 2.85 -0.17 (21) 195410 56.50 1.15 2.28 2.08 -0.20 (22) 195411 41.70 5.15 0.60 0.53 -0.07 (23) 195412 21.00 0.34 0.00 0.00 0.00 (24) 195501 20.60 1.53 0.00 0.00 0.00 (25) 195502 28.50 1.44 0.00 0.00 0.00 (26) 195503 39.90 1.46 0.62 0.51 -0.11 (27) 195504 46.40 1.24 1.39 1.22 -0.17 (28) 195505 54.80 4.13 2.70 2.50 -0.20 (29) 195506 65.40 8.65 4.27 4.14 -0.13 (30) 195507 78.60 4.43 6.66 6.51 -0.15 (31) 195508 72.90 1.64 5.25 5.12 -0.13 (32) 195509 64.50 3.89 3.46 3.30 -0.16 (33) 195510 50.70 3.55 1.60 1.46 -0.14 (34) 195511 33.90 2.89 0.07 0.06 -0.01 (35) 195512 21.60 1.31 0.00 0.00 0.00Example 3
Read an ascii file containing 117 years of monthly temperatures at Boulder, Colorado (source: http://www.esrl.noaa.gov/psd/boulder/). Compute potential evapotranspiration (mm/month).
dirt = "./" filt = "Boulder.temp.1897-2013.txt" ncol = 14 ; due to the structure of the data, it is better to read as 2D nrow = numAsciiRow(dirt+filt) ; contributed.ncl tmp2 = asciiread(dirt+filt,(/nrow,ncol/), "float") tmp = ndtooned(tmp2(:,1:ncol-2)) ; make one dimensional tmp@units = "degF" tmp@long_name = "Boulder Temperature" tmp = (5.0/9.0)*(tmp-32) ; Convert F=>C tmp@units = "degC" printVarSummary(tmp) print("tmp: min="+min(tmp)+" max="+max(tmp)) TMP = ndtooned( tmp ) ; create a time series ; latBld = 40.0 ; Boulder is at 40N latitude pet = thornthwaite(TMP, latBld, False, 0) pet@long_name = "PET" pet@units = "mm/month" ; ie, monthly total printVarSummary(pet) print("pet: min="+min(tmp)+" max="+max(tmp)) write_matrix (onedtond(pet,(/nyrs,12/)), "12f7.1", False) ; requires 2DThe (edited) outputVariable: pet Type: float Total Size: 11232 bytes 1404 values Number of Dimensions: 1 Dimensions and sizes: [1404] Coordinates: Number Of Attributes: 3 units : mm/month long_name : PET: Boulder, CO _FillValue : -9999 (0) tmp: min=-8.66667 max=25.7222 --- Variable: pet Type: float Total Size: 5616 bytes 1404 values Number of Dimensions: 1 Dimensions and sizes: [1404] Coordinates: Number Of Attributes: 3 units : mm/month long_name : PET _FillValue : 9.96921e+36 (0) pet: min=0 max=162.275 --- J F M A M J J A S O N D 0.0 0.2 7.4 38.6 86.6 103.7 129.9 115.6 94.0 44.4 16.0 0.0 0.0 10.1 6.0 42.2 56.5 111.6 136.5 131.1 78.7 35.9 6.3 0.0 0.0 0.0 3.0 41.9 74.7 108.8 123.3 121.0 91.3 39.9 28.4 4.9 7.1 0.0 23.1 32.2 85.1 123.5 130.3 128.0 77.1 58.2 21.7 10.2 6.0 0.0 15.3 36.0 83.9 109.2 149.1 124.5 76.4 48.2 28.7 3.7 0.1 8.7 13.2 41.0 82.8 112.4 124.1 126.4 74.5 52.0 16.4 6.0 [SNIP] 0.0 5.5 17.2 38.0 75.1 109.6 148.2 117.9 74.8 28.3 24.7 0.0 9.2 11.8 25.9 36.6 83.2 97.5 125.8 117.9 81.9 24.7 20.1 0.0 1.0 0.0 21.9 41.0 63.5 112.8 137.7 128.8 93.3 52.0 12.1 7.2 1.1 0.0 28.3 41.3 62.8 115.5 141.9 139.4 82.8 46.3 17.0 0.1 10.5 0.2 43.6 58.0 86.6 142.3 147.4 131.9 91.3 40.7 24.9 1.8 1.0 0.1 16.5 26.6 77.3 124.7 136.5 128.0 88.4 32.9 18.8 0.0Example 4
Read the Wichita data used by the R-SPEI package. Calculate the Thornthwaite estimates via NCL's thornthwaite function and compare with those on the file. NOTE: The values are basically identical.
dirw = "./" filw = "Wichita.Rspei.198001-201112.nc" fw = addfile(dirw+filw, "r") tavg = fw->TAVG printVarSummary(tavg) print("tavg: min="+min(tavg)+" max="+max(tavg)) lat = fw@latitude ; Wichita: 37.65N pet = thornthwaite(tavg, lat, False, 0) pet@long_name = "PET: Wichita" pet@units = "mm/month" ; ie, monthly total pet@tag = "NCL: thormthwaite" copy_VarCoords(tavg, pet) printVarSummary(pet) print("pet: min="+min(pet)+" max="+max(pet)) ntim = dimsizes(tavg) nyrs = ntim/12 write_matrix (onedtond(pet,(/nyrs,12/)), "12f7.1", False) ; Read SPEI-R values from R-SPEI file for comparison. ; Note the last 2 months contain missing values. pet_spei = fw->PET_TW pet_diff = max( abs(pet - pet_spei) ) print("maximum difference: pet_diff="+pet_diff)The output is:
Variable: tavg Type: float Total Size: 1536 bytes 384 values Number of Dimensions: 1 Dimensions and sizes: [time | 384] Coordinates: time: [198001..201112] Number Of Attributes: 3 long_name : average temperature units : degC _FillValue : -999 (0) tavg: min=-8.71 max=32.46 Variable: pet Type: float Total Size: 1536 bytes 384 values Number of Dimensions: 1 Dimensions and sizes: [time | 384] Coordinates: time: [198001..201112] Number Of Attributes: 4 tag : NCL: thormthwaite units : mm/month long_name : PET: Wichita _FillValue : -999 (0) pet: min=0 max=226.929 J F M A M J J A S O N D 0.0 0.0 10.9 44.2 84.4 163.2 226.9 184.8 117.6 51.7 18.3 3.1 0.8 6.9 23.7 76.9 80.7 152.6 186.5 147.6 104.2 43.2 18.3 0.2 [SNIP] 0.0 0.1 22.6 66.6 94.2 167.7 181.6 175.5 113.6 62.8 17.4 0.9 0.0 0.1 24.6 58.7 97.2 171.7 220.2 183.3 93.7 81.3 -999.0 -999.0 (0) maximum difference: pet_diff=1.52588e-05