全国580+ 站点，20分钟全搞定。

详情邮件联系：liyj AT msu.edu

部分代码如下：

mutate(Jd      = dayid,

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

                                                        # h: elevation above sea level [m]

           lat_rad = deg2rad(lat), 

           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),

           ws      = acos(-tan(lat_rad)*tan(delta_low)),

           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      = (24*60*0.082*dr/3.1415926)*(ws*sin(lat_rad)*sin(delta_low)+cos(lat_rad)*cos(delta_low)*sin(ws)),

                           # 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]

    )