# 全国ETo (FAO Penman-Monteith) 1960-2017逐日计算 R代码

mutate(Jd      = dayid,

P       = 101.3*(((293-0.0065*h)/293)^5.26), # P: atmospheric pressure [kPa]

# h: elevation above sea level [m]

cp      = 1.013*10^(-3),     # cp: specific heat at constant pressure, 1.013 10-3 [MJ kg-1 °C-1]

e       = 0.622,             # e: ratio molecular weight of water vapour/dry air = 0.622

lam     = 2.45,              # lam: latent heat of vaporization, 2.45 [MJ kg-1],

r       = cp*P/e/lam,        # psychrometric constant [kPa °C-1]

Tmean   = (tmax+tmin)/2,     # mean daily air temperature at 2 m height [°C]

e0T     = 0.6108*exp(17.27*Tmean/(Tmean+237.3)), # e°(T) saturation vapour pressure at the air temperature T [kPa]

e0Tmax  = 0.6108*exp(17.27*tmax/(tmax+237.3)),

e0Tmin  = 0.6108*exp(17.27*tmin/(tmin+237.3)),

es      = (e0Tmax+e0Tmin)/2, # saturation vapour pressure [kPa]

ea      = (rh/100)*e0T,      # actual vapour pressure [kPa]

# es - ea saturation vapour pressure deficit [kPa]

delta_low = 0.409*sin(2*3.1415927*Jd/365-1.39),

n       = tsun,              # n actual duration of sunshine [hour]

N       = (24/3.1415927)*ws, # N maximum possible duration of sunshine or daylight hours [hour]

dr      = 1+0.033*cos(2*3.1415927*Jd/365),

# Ra extraterrestrial radiation [MJ m-2 day-1]

a       = 0.25, # empirical coefficient 0.12~0.29

b       = 0.50, # empirical coefficient 0.45~0.73

Rs      = (a + b*n/N)*Ra, # Rs solar or shortwave radiation [MJ m-2 day-1]

Rns     = (1-0.23)*Rs,

Rso     = (0.75+2*(10)^(-5)*h)*Ra,

Tkm     = (tmax+273.2),

Tkn     = (tmin+273.2),

Rnl     = 4.903*(10^(-9))*(0.1+0.9*n/N)*(0.34-0.14*(ea^(0.5)))*((Tkm^4+Tkn^4)/2),

Rn      = Rns - Rnl,                 # net radiation at the crop surface [MJ m-2 day-1]

Delta   = 4098*es/((Tmean+237.3)^2), # slope vapour pressure curve [kPa °C-1]

uz      = wind * 0.72,       # uz measured wind speed at z m above ground surface [m s-1]

# should multiply 0.72 (empirical coefficient)to correct

z       = 10,                # z height of measurement above ground surface [m]

# measurement at 10 m (Changming Liu)

u2      = 4.87*uz/log(67.8*z-5.42), # u2 wind speed at 2 m height [m s-1]

G       = 0,      # soil heat flux density [MJ m-2 day-1]

# For day and ten-day period:    Gday = 0

# For hourly or shorter periods: Ghr  = 0.1*Rn (daylight periods); Ghr = 0.5*Rn (during nighttime period)

ET0     = (0.408*Delta*(Rn-G) + 900*r*u2*(es-ea)/(Tmean+273)) / (Delta+r*(1+0.34*u2)) # reference evapotranspiration [mm day-1]

)

http://blog.sciencenet.cn/blog-2073314-1138140.html

## 全部精选博文导读

GMT+8, 2020-10-24 19:49