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

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

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