Nothing
# R version 3.2.3 (2015-12-10) -- "Wooden Christmas-Tree"
#
# evap.R
#
# VERSION: 1.0-r2
# LAST UPDATED: 2016-09-11
#
# ~~~~~~~~
# license:
# ~~~~~~~~
# Copyright (C) 2016 Prentice Lab
#
# This file is part of the SPLASH model.
#
# SPLASH is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 2.1 of the License, or
# (at your option) any later version.
#
# SPLASH is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with SPLASH. If not, see <http://www.gnu.org/licenses/>.
#
# ~~~~~~~~~
# citation:
# ~~~~~~~~~
# T. W. Davis, I. C. Prentice, B. D. Stocker, R. J. Whitley, H. Wang, B. J.
# Evans, A. V. Gallego-Sala, M. T. Sykes, and W. Cramer, Simple process-led
# algorithms for simulating habitats (SPLASH): Robust indices of radiation
# evapo-transpiration and plant-available moisture, Geoscientific Model
# Development, 2016 (in progress)
#
# ~~~~~~~~~~~~
# description:
# ~~~~~~~~~~~~
# This script contains functions to calculate daily radiation, condensation,
# and evapotranspiration, i.e.:
# berger_tls(double n, double N)
# density_h2o(double tc, double pa)
# dcos(double d)
# dsin(double d)
# elv2pres(double z)
# enthalpy_vap(double tc)
# evap(double lat, double n, double elv=0, double y=0, double sf=1,
# double tc=23.0, double sw=1.0)
# julian_day(double y, double m, double i)
# psychro(double tc, double pa)
# sat_slope(double tc)
# specific_heat(double tc)
#
# ~~~~~~~~~~
# changelog:
# ~~~~~~~~~~
# - fixed Cooper's and Spencer's declination angle equations [14.11.25]
# - replaced simplified_kepler with full_kepler method [14.11.25]
# - added berger_tls function [15.01.13]
# - updated evap function (similar to stash.py EVAP class) [15.01.13]
# - added missing variables to evap list [16.02.17]
# - updated documentation [16.05.27]
# - addressed specific heat limitation [16.09.11]
#
#### IMPORT SOURCES ##########################################################
# source("const.R")
# source("solar.R")
#### DEFINE FUNCTIONS ########################################################
# ************************************************************************
# Name: density_h2o
# Inputs: - double (tc), air temperature, degrees C
# - double (pa), atm pressure, Pa
# Returns: double, kg/m^3
# Features: This function calculates the temperature and pressure
# dependent density of pure water
# * Ref: Chen, C.T., R.A. Fine, and F.J. Millero (1977), The equation
# of state of pure water determined from sound speeds, The
# Journal of Chemical Physics 66, 2142;
# \doi{10.1063/1.434179}
# ************************************************************************
#' Calculate density of water at 1 atm, g/cm^3
#'
#' This function calculates the temperature and pressure dependent density of
#' pure water.
#'
#' @param tc double, air temperature, degrees C.
#' @param pa double, atm pressure, Pa.
#'
#' @return double, kg/m^3.
#' @keywords internal
#'
#' @references
#' Chen, C.T., Fine, R.A. and Millero, F.J., 1977. The equation of state of
#' pure water determined from sound speeds. The Journal of Chemical Physics,
#' 66(5), pp.2142-2144. \doi{10.1063/1.434179}
density_h2o <- function(tc, pa) {
# Calculate density of water at 1 atm, g/cm^3
po <- 0.99983952 +
(6.788260e-5) * tc +
-(9.08659e-6) * tc * tc +
(1.022130e-7) * tc * tc * tc +
-(1.35439e-9) * tc * tc * tc * tc +
(1.471150e-11) * tc * tc * tc * tc * tc +
-(1.11663e-13) * tc * tc * tc * tc * tc * tc +
(5.044070e-16) * tc * tc * tc * tc * tc * tc * tc +
-(1.00659e-18) * tc * tc * tc * tc * tc * tc * tc * tc
# Calculate the bulk modulus of water at 1 atm, atm
ko <- 19652.17 +
148.1830 * tc +
-2.29995 * tc * tc +
0.01281 * tc * tc * tc +
-(4.91564e-5) * tc * tc * tc * tc +
(1.035530e-7) * tc * tc * tc * tc * tc
# Calculate temperature-dependend coefficients
ca <- 3.26138 +
(5.223e-4) * tc +
(1.324e-4) * tc * tc +
-(7.655e-7) * tc * tc * tc +
(8.584e-10) * tc * tc * tc * tc
cb <- (7.2061e-5) +
-(5.8948e-6) * tc +
(8.69900e-8) * tc * tc +
-(1.0100e-9) * tc * tc * tc +
(4.3220e-12) * tc * tc * tc * tc
# Convert pressure to bar (1 bar = 100000 Pa)
pbar <- (1e-5) * pa
pw <-
(1e3) * po * (ko + ca * pbar + cb * pbar ^ 2) / (ko + ca * pbar + cb *
pbar ^ 2 - pbar)
return(pw)
}
# ************************************************************************
# Name: elv2pres
# Inputs: double (z), meters
# Returns: double, Pa
# Features: Calculates atmospheric pressure for a given elevation
# Depends: - kPo ............ base pressure, Pa
# - kTo ............ base temperature, C
# - kL ............. temperature lapse rate, K/m
# - kR ............. universal gas constant, J/mol/K
# - kMa ............ molecular weight of dry air, kg/mol
# - kG ............. gravity, m/s^2
# Ref: Allen et al. (1998)
# ************************************************************************
#' Elevation to pressure
#'
#' Calculates atmospheric pressure for a given elevation.
#'
#' @param z Elevation, m.
#' @param kG Gravitational acceleration, m/s^2.
#' Default: \eqn{9.80665} (Allen, 1973)
#' @param kL Adiabatic lapse rate, K/m.
#' Default: \eqn{0.0065} (Allen, 1973)
#' @param kMa Molecular weight of dry air, kg/mol.
#' Default: \eqn{0.028963} (Tsilingiris, 2008)
#' @param kPo Standard atmosphere, Pa.
#' Default: \eqn{101325} (Allen, 1973)
#' @param kR Universal gas constant, J/mol/K.
#' Default: \eqn{8.31447} (Moldover et al., 1988)
#' @param kTo Base temperature, K.
#' Default: \eqn{288.15} (Berberan-Santos et al., 1997)
#'
#' @return Atmospheric pressure for the given elevation, Pa.
#' @keywords internal
#'
#' @references
#' Allen, R.G., Pereira, L.S., Raes, D. and Smith, M., 1998. Crop
#' evapotranspiration-Guidelines for computing crop water requirements-FAO
#' Irrigation and drainage paper 56. Food and Agriculture Organization of the
#' United Nations, Rome, 300(9), p.D05109. Available:
#' http://www.fao.org/docrep/x0490e/x0490e07.htm
#'
#' Berberan-Santos, M.N., Bodunov, E.N. and Pogliani, L., 1997. On the
#' barometric formula. American Journal of Physics, 65(5), pp.404-412.
#' \doi{10.1119/1.18555}
#'
#' Moldover, M.R., Trusler, J.M., Edwards, T.J., Mehl, J.B. and Davis, R.S.,
#' 1988. Measurement of the universal gas constant R using a spherical acoustic
#' resonator. Physical review letters, 60(4), p.249.
#' \doi{10.1103/PhysRevLett.60.249}
#'
#' Tsilingiris, P.T., 2008. Thermophysical and transport properties of humid air
#' at temperature range between 0 and 100 C. Energy Conversion and Management,
#' 49(5), pp.1098-1110. \doi{10.1016/j.enconman.2007.09.015}
elv2pres <- function(z,
kG = 9.80665,
kL = 0.0065,
kMa = 0.028963,
kPo = 101325,
kR = 8.31447,
kTo = 288.15) {
kPo * (1 - kL * z / kTo) ^ (kG * kMa / (kR * kL))
}
# ************************************************************************
# Name: enthalpy_vap
# Inputs: double (tc), air temperature, degrees C
# Returns: double, J/kg
# Features: This function calculates the temperature-dependent enthalpy
# of vaporization (latent heat of vaporization)
# Ref: Eq. 8, Henderson-Sellers (1984), A new formula for latent heat
# of vaporization of water as a function of temperature, Quarterly
# Journal of the Royal Meteorological Society, vol. 110, pp. 1186--
# 1190.
# ************************************************************************
#' Calculate enthalpy of vaporization
#'
#' This function calculates the temperature-dependent enthalpy of vaporization
#' (latent heat of vaporization).
#'
#' @inheritParams density_h2o
#'
#' @return double, J/kg.
#' @keywords internal
#'
#' @references
#' Eq. 8, Henderson‐Sellers, B., 1984. A new formula for latent heat of
#' vaporization of water as a function of temperature. Quarterly Journal of the
#' Royal Meteorological Society, 110(466), pp.1186-1190.
#' \doi{10.1002/qj.49711046626}
enthalpy_vap <- function(tc) {
1.91846e6 * ((tc + 273.15) / (tc + 273.15 - 33.91)) ^ 2
}
# ************************************************************************
# Name: psychro
# Inputs: - double (tc), air temperature, degrees C
# - double (pa), atm pressure, Pa
# Returns: double, Pa/K
# Features: This function calculates the temperature and pressure
# dependent psychrometric constant
# Depends: - enthalpy_vap
# - specific_heat
# Ref: Allen, R.G., L.S. Pereira, D. Raes, M. Smith (1998),
# 'Meteorological data,' Crop evapotranspiration - Guidelines
# for computing crop water requirements - FAO Irrigation and
# drainage paper 56, Food and Agriculture Organization of the
# United Nations, Available:
# http://www.fao.org/docrep/x0490e/x0490e07.htm
# ************************************************************************
#' Calculate psychrometric constant
#'
#' This function calculates the temperature and pressure dependent
#' psychrometric constant.
#'
#' @inheritParams density_h2o
#' @param kMa double, molecular weight of dry air, kg/mol.
#' Default: \eqn{0.028963} (Tsilingiris, 2008)
#' @param kMv double, molecular weight of water vapor, kg/mol.
#' Default: \eqn{0.01802} (Tsilingiris, 2008)
#' @return double, Pa/K.
#' @keywords internal
#'
#' @references
#' Allen, R.G., Pereira, L.S., Raes, D. and Smith, M., 1998. Crop
#' evapotranspiration-Guidelines for computing crop water requirements-FAO
#' Irrigation and drainage paper 56. Food and Agriculture Organization of the
#' United Nations, Rome, 300(9), p.D05109. Available:
#' http://www.fao.org/docrep/x0490e/x0490e07.htm
psychro <- function(tc, pa, kMa = 0.028963, kMv = 0.01802) {
# Calculate the specific heat capacity of water, J/kg/K
cp <- specific_heat(tc)
# Calculate latent heat of vaporization, J/kg
lv <- enthalpy_vap(tc)
# Calculate psychrometric constant, Pa/K
return(cp * kMa * pa / (kMv * lv))
}
# ************************************************************************
# Name: specific_heat
# Inputs: double (tc), air temperature, degrees C
# Returns: double, specific heat of moist air, J/kg/K
# Features: This function calculates the spefic heat of moist air
# Ref: Tsilingris (2008), Thermophysical and transport properties of
# humid air at temperature range between 0 and 100 °C, Energy
# Conversion and Management, vol. 49, pp. 1098--1110.
# ************************************************************************
#' Calculate specific heat
#'
#' This function calculates the specific heat of moist air.
#'
#' @inheritParams density_h2o
#'
#' @return double, specific heat of moist air, J/kg/K.
#' @keywords internal
#'
#' @references
#' Tsilingiris, P.T., 2008. Thermophysical and transport properties of humid air
#' at temperature range between 0 and 100 C. Energy Conversion and Management,
#' 49(5), pp.1098-1110. \doi{10.1016/j.enconman.2007.09.015}
specific_heat <- function(tc) {
if (tc < 0) {
tc <- 0
} else if (tc > 100) {
tc <- 100
}
cp <- 1.0045714270 +
(2.050632750e-3) * tc -
(1.631537093e-4) * tc * tc +
(6.212300300e-6) * tc * tc * tc -
(8.830478888e-8) * tc * tc * tc * tc +
(5.071307038e-10) * tc * tc * tc * tc * tc
cp <- (1e3) * cp
return(cp)
}
# ************************************************************************
# Name: sat_slope
# Inputs: double (tc), degrees C
# Returns: double, Pa/K
# Features: This function calculates the temperature-dependent slope of
# the saturation pressure temperature curve using the
# methodology presented in the eMast energy.cpp script
# Ref: - Eq. 6, Prentice et al. (1993);
# - Eq. 13, Allen et al. (1998)
# ************************************************************************
#' Calculate the temperature-dependent slope
#'
#' This function calculates the temperature-dependent slope of the saturation
#' pressure temperature curve using the methodology presented in the eMast
#' energy.cpp script.
#' @inheritParams density_h2o
#'
#' @return double, Pa/K.
#' @keywords internal
#'
#' @references
#' Allen, R.G., Pereira, L.S., Raes, D. and Smith, M., 1998. Crop
#' evapotranspiration-Guidelines for computing crop water requirements-FAO
#' Irrigation and drainage paper 56. Food and Agriculture Organization of the
#' United Nations, Rome, 300(9), p.D05109. Available:
#' http://www.fao.org/docrep/x0490e/x0490e07.htm
#'
#' Prentice, I.C., Sykes, M.T. and Cramer, W., 1993. A simulation model for the
#' transient effects of climate change on forest landscapes. Ecological
#' modelling, 65(1-2), pp.51-70. \doi{10.1016/0304-3800(93)90126-D}
sat_slope <- function(tc) {
(17.269) * (237.3) * (610.78) *
exp(17.269 * tc / (237.3 + tc)) / (237.3 + tc) ^ 2
}
# ************************************************************************
# Name: calc_daily_evap
# Inputs: - double, latitude, degrees (lat)
# - double, day of year (n)
# - double, elevation (elv) *optional
# - double, year (y) *optional
# - double, fraction of sunshine hours (sf) *optional
# - double, mean daily air temperature, deg C (tc) *optional
# - double, evaporative supply rate, mm/hr (sw) *optional
# Returns: list object (evap)
# $nu_deg ............ true anomaly, degrees
# $lambda_deg ........ true longitude, degrees
# $dr ................ distance factor, unitless
# $delta_deg ......... declination angle, degrees
# $hs_deg ............ sunset angle, degrees
# $ra_j.m2 ........... daily extraterrestrial radiation, J/m^2
# $tau ............... atmospheric transmittivity, unitless
# $ppfd_mol.m2 ....... daily photosyn photon flux density, mol/m^2
# $hn_deg ............ net radiation hour angle, degrees
# $rn_j.m2 ........... daily net radiation, J/m^2
# $rnn_j.m2 .......... daily nighttime net radiation, J/m^2
# $econ_m3.j ......... water to energy conversion, m^3/J
# $cond_mm ........... daily condensation, mm
# $eet_mm ............ daily equilibrium evapotranspiration, mm
# $pet_mm ............ daily potential evapotranspiration, mm
# $hi_deg ............ intersection hour angle, degrees
# $aet_mm ............ daily actual evapotranspiration, mm
# Features: This function calculates daily radiation, condensation, and
# evaporation fluxes.
# Depends: - kw ............. entrainment factor for PET
# - calc_daily_solar daily radiation fluxes
# - dcos() ......... cos(x*pi/180), where x is in degrees
# - dsin() ......... sin(x*pi/180), where x is in degrees
# - density_h2o() .. density of water
# - elv2pres() ..... elevation dependent atm. pressure
# - enthalpy_vap() . latent heat of vaporization
# - psychro() ...... psychrometric constant
# - sat_slope() .... slope of sat. pressure temp curve
# ************************************************************************
#' Calculate daily evaporation fluxes
#'
#' This function calculates daily radiation, condensation, and evaporation
#' fluxes.
#'
#' @param lat double, decimal degrees.
#' @param n double, day of year.
#' @param elv double, elevation, m A.S.L.
#' Default: \eqn{0}.
#' @param y double, year.
#' Default: \eqn{0}.
#' @param sf double, fraction of sunshine hours.
#' Default: \eqn{1}.
#' @param tc double, mean daily air temperature, degrees C.
#' Default: \eqn{23.0}.
#' @param sw double, evaporative supply rate, mm/hr.
#' Default: \eqn{1.0}.
#' @param ke double, eccentricity of earth's orbit.
#' Default: \eqn{0.01670}, 2000CE (Berger, 1978).
#' @param keps double, obliquity of earth's elliptic.
#' Default: \eqn{23.44}, 2000CE (Berger, 1978).
#' @param komega double, lon. of perihelion, degrees
#' Default: \eqn{283}, 2000CE (Berger, 1978).
#' @param kw double, PET entrainment, \eqn{(1 + kw) * EET}
#' Default: \eqn{0.26} (Priestley-Taylor, 1972)
#'
#' @return Returns a \code{list} object with the following variables:
#' \itemize{
#' \item nu_deg ............ true anomaly, degrees
#' \item lambda_deg ........ true longitude, degrees
#' \item dr ................ distance factor, unitless
#' \item delta_deg ......... declination angle, degrees
#' \item hs_deg ............ sunset angle, degrees
#' \item ra_j.m2 ........... daily extraterrestrial radiation, J/m^2
#' \item tau ............... atmospheric transmittivity, unitless
#' \item ppfd_mol.m2 ....... daily photosyn photon flux density, mol/m^2
#' \item hn_deg ............ net radiation hour angle, degrees
#' \item rn_j.m2 ........... daily net radiation, J/m^2
#' \item rnn_j.m2 .......... daily nighttime net radiation, J/m^2
#' \item econ_m3.j ......... water to energy conversion, m^3/J
#' \item cond_mm ........... daily condensation, mm
#' \item eet_mm ............ daily equilibrium evapotranspiration, mm
#' \item pet_mm ............ daily potential evapotranspiration, mm
#' \item hi_deg ............ intersection hour angle, degrees
#' \item aet_mm ............ daily actual evapotranspiration, mm
#' }
#' @export
#'
#' @references
#' Berger, A.L., 1978. Long-term variations of daily insolation and Quaternary
#' climatic changes. Journal of Atmospheric Sciences, 35(12), pp.2362-2367.
#' \doi{10.1175/1520-0469(1978)035<2362:LTVODI>2.0.CO;2}
#'
#' Priestley, C.H.B. and Taylor, R.J., 1972. On the assessment of surface heat
#' flux and evaporation using large-scale parameters. Monthly weather review,
#' 100(2), pp.81-92. \doi{10.1175/1520-0493(1972)100<0081:OTAOSH>2.3.CO;2}
#'
#' @examples
#' evap <- splash::calc_daily_evap(lat = 37.7,
#' n = 172,
#' elv = 142,
#' y = 2000,
#' sf = 1,
#' tc = 23.0,
#' sw = 0.9)
#' cat(sprintf("Evaporation values:\n"))
#' cat(sprintf(" s: %0.6f Pa/K\n", evap$s_pa.k))
#' cat(sprintf(" Lv: %0.6f MJ/kg\n", (1e-6) * evap$lv_j.kg))
#' cat(sprintf(" Patm: %0.6f bar\n", (1e-5) * evap$patm_pa))
#' cat(sprintf(" pw: %0.6f kg/m^3\n", evap$pw_kg.m3))
#' cat(sprintf(" gamma: %0.6f Pa/K\n", evap$gam_pa.k))
#' cat(sprintf(" Econ: %0.6f mm^3/J\n", (1e9) * evap$econ_m3.j))
#' cat(sprintf(" Cn: %0.6f mm\n", evap$cond_mm))
#' cat(sprintf(" rx: %0.6f\n", evap$rx))
#' cat(sprintf(" hi: %0.6f degrees\n", evap$hi_deg))
#' cat(sprintf(" EET: %0.6f mm\n", evap$eet_mm))
#' cat(sprintf(" PET: %0.6f mm\n", evap$pet_mm))
#' cat(sprintf(" AET: %0.6f mm\n", evap$aet_mm))
calc_daily_evap <-function(lat,
n,
elv = 0,
y = 0,
sf = 1,
tc = 23.0,
sw = 1.0,
ke = 0.01670,
keps = 23.44,
komega = 283,
kw = 0.26) {
# Local bindings
pir <- pi / 180
# ~~~~~~~~~~~~~~~~~~~~~~~~ FUNCTION WARNINGS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
if (lat > 90 || lat < -90) {
stop("Warning: Latitude outside range of validity (-90 to 90)!")
}
if (n < 1 || n > 366) {
stop("Warning: Day outside range of validity (1 to 366)!")
}
# ~~~~~~~~~~~~~~~~~~~~~~~ FUNCTION VARIABLES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
evap <- list()
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 1. Calculate radiation fluxes
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
solar <- calc_daily_solar(lat = lat,
n = n,
elv = elv,
y = y,
sf = sf,
tc = tc,
ke = ke,
keps = keps,
komega = komega)
ru <- solar$ru;
rv <- solar$rv;
rw <- solar$rw;
rnl <- solar$rnl_w.m2;
hn <- solar$hn_deg;
rn_d <- solar$rn_j.m2;
rnn_d <- solar$rnn_j.m2;
evap$ra_j.m2 <- solar$ra_j.m2
evap$rn_j.m2 <- solar$rn_j.m2
evap$ppfd_mol.m2 <- solar$ppfd_mol.m2
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 2. Calculate water-to-energy conversion (econ), m^3/J
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Atmospheric pressure, Pa
patm <- elv2pres(elv)
evap$patm_pa <- patm
# Slope of saturation vap press temp curve, Pa/K
s <- sat_slope(tc)
evap$s_pa.k <- s
# Enthalpy of vaporization, J/kg
lv <- enthalpy_vap(tc)
evap$lv_j.kg <- lv
# Density of water, kg/m^3
pw <- density_h2o(tc, patm)
evap$pw_kg.m3 <- pw
# Psychrometric constant, Pa/K
gam <- psychro(tc, patm)
evap$gam_pa.k <- gam
econ <- s / (lv * pw * (s + gam))
evap$econ_m3.j <- econ
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 3. Calculate daily condensation (cn), mm
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cn <- (1e3) * econ * abs(rnn_d)
evap$cond_mm <- cn
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 4. Estimate daily equilibrium evapotranspiration (eet_d), mm
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
eet_d <- (1e3) * econ * rn_d
evap$eet_mm <- eet_d
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 5. Estimate daily potential evapotranspiration (pet_d), mm
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pet_d <- (1 + kw) * eet_d
evap$pet_mm <- pet_d
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 6. Calculate variable substitute (rx), (mm/hr)/(W/m^2)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rx <- (3.6e6) * (1 + kw) * econ
evap$rx <- rx
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 7. Calculate the intersection hour angle (hi), degrees
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cos_hi <- sw / (rw * rv * rx) + rnl / (rw * rv) - ru / rv
# print( paste( "in evap: sw =", sw))
# print( paste( "in evap: ru =", ru))
# print( paste( "in evap: rv =", rv))
# print( paste( "in evap: rw =", rw))
# print( paste( "in evap: rx =", rx))
# print( paste( "in evap: rnl =", rnl))
# print( paste( "in evap: cos_hi =", cos_hi))
if (cos_hi >= 1.0) {
hi <- 0.0 # supply exceeds demand
} else if (cos_hi <= -1.0) {
hi <- 180.0 # supply limits demand everywhere
} else {
hi <- acos(cos_hi)
hi <- hi / pir
}
evap$hi_deg <- hi
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 8. Estimate daily actual evapotranspiration (aet_d), mm
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
aet_d <- (24 / pi) * (
sw * hi * pir +
rx * rw * rv * (dsin(hn) - dsin(hi)) +
(rx * rw * ru - rx * rnl) * (hn - hi) * pir
)
evap$aet_mm <- aet_d
# ~~~~~~~~~~~~~~~~~~~~~~~~~~ RETURN VALUES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #
return(evap)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.