NCL Home > Documentation > Functions > Crop, Meteorology

# 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

t

Monthly 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).

lat

Latitude(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

opt

Not used. Set to False.

dim

The 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

```

## 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.00

Example 3
Read an ascii file containing 117 years of monthly temperatures
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

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 2D

The (edited) output

Variable: 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.0

Example 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"

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

```