lyjcas的个人博客分享 http://blog.sciencenet.cn/u/lyjcas

博文

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

已有 1843 次阅读 2018-10-1 08:50 |个人分类:data and tools|系统分类:科研笔记

 


全国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]

    )

  




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

上一篇:Ph.D. career
下一篇:[招募博士研究生合作研究] 大流域社会经济与自然系统耦合项目

1 梁守真

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-11-24 18:26

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部