successful test at single site: correct forcing, p-model runs

DESCRIPTION

@@ -13,6 +13,7 @@ Description: A wrapper for global rsofun runs
 License: AGPL-3
 Encoding: UTF-8
 Imports: 
+    zoo,
     dplyr,
     tidyr,
     magrittr,

R/calc_patm.R

@@ -0,0 +1,47 @@
+#' Calculates atmospheric pressure
+#'
+#' Calculates atmospheric pressure as a function of elevation, by default assuming
+#' standard atmosphere (101325 Pa at sea level)
+#'
+#' @param elv Elevation above sea-level (m.a.s.l.)
+#' @param patm0 (Optional) Atmospheric pressure at sea level (Pa), defaults to 101325 Pa.
+#'
+#' @details The elevation-dependence of atmospheric pressure is computed by
+#' assuming a linear decrease in temperature with elevation and a mean
+#' adiabatic lapse rate (Berberan-Santos et al., 1997):
+#' \deqn{
+#'    p(z) = p0 ( 1 - Lz / TK0) ^ ( g M / (RL) )
+#' }
+#' where \eqn{z} is the elevation above mean sea level (m, argument \code{elv}),
+#' \eqn{g} is the gravity constant (9.80665 m s-2), \eqn{p0} is the atmospheric
+#' pressure at 0 m a.s.l. (argument \code{patm0}, defaults to 101325 Pa),
+#' \eqn{L} is the mean adiabatic lapse rate (0.0065 K m-2),
+#' \eqn{M} is the molecular weight for dry air (0.028963 kg mol-1),
+#' \eqn{R} is the universal gas constant (8.3145 J mol-1 K-1), and \eqn{TK0}
+#' is the standard temperature (298.15 K, corresponds to 25 deg C).
+#'
+#' @return A numeric value for \eqn{p}
+#'
+#' @examples print("Standard atmospheric pressure, in Pa, corrected for 1000 m.a.s.l.:")
+#' print(calc_patm(1000))
+#'
+#' @references  Allen, R. G., Pereira, L. S., Raes, D., Smith, M.:
+#'              FAO Irrigation and Drainage Paper No. 56, Food and
+#'              Agriculture Organization of the United Nations, 1998
+#'
+#' @export
+#'
+calc_patm <- function( elv, patm0 = 101325 ){
+
+  # Define constants:
+  kTo <- 298.15    # base temperature, K (Prentice, unpublished)
+  kL  <- 0.0065    # adiabiatic temperature lapse rate, K/m (Allen, 1973)
+  kG  <- 9.80665   # gravitational acceleration, m/s^2 (Allen, 1973)
+  kR  <- 8.3145    # universal gas constant, J/mol/K (Allen, 1973)
+  kMa <- 0.028963  # molecular weight of dry air, kg/mol (Tsilingiris, 2008)
+
+  # Convert elevation to pressure, Pa:
+  patm <- patm0*(1.0 - kL*elv/kTo)^(kG*kMa/(kR*kL))
+
+  return(patm)
+}

R/calc_vp.R

@@ -4,49 +4,49 @@
 #'
 #' @param qair Air specific humidity (g g-1)
 #' @param patm Atmospehric pressure (Pa)
-#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing 
+#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing
 #' for calculating it based on standard sea level pressure)
 #'
 #' @return vapor pressure (Pa)
 #' @export
-#' 
+#'
 calc_vp <- function(
   qair,
   patm = NA,
   elv = NA
 ) {
-  
-  # Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure 
-  #           Calculation Methods, in Evaporation and Evapotranspiration: 
+
+  # Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure
+  #           Calculation Methods, in Evaporation and Evapotranspiration:
   #           Measurements and Estimations, Springer, London.
   #             vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea
   #             where:
   #                 tc = average daily air temperature, deg C
   #                 eact  = actual vapor pressure, Pa
-
 
-  # calculate atmopheric pressure (Pa) assuming standard 
+
+  # calculate atmopheric pressure (Pa) assuming standard
   # conditions at sea level (elv=0)
   # if (is.na(elv) && is.na(patm)){
-  #   
+  #
   #   warning(
   #     "
-  #     calc_vp(): Either patm or elv must be provided 
+  #     calc_vp(): Either patm or elv must be provided
   #     if eact is not given.
   #     ")
   #   vp <- NA
-  #   
+  #
   # } else {
 
     patm <- ifelse(is.na(patm),
                    calc_patm(elv),
                    patm)
-    
+
     # Calculate VP.
     vp <- calc_vp_inst(qair = qair, patm = patm)
   # }
   return( vp )
-  
+
 }
 
 #' Instantaneous vapour pressure
@@ -63,26 +63,26 @@ calc_vp_inst <- function(
   qair,
   patm
 ){
-  
-  # Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure 
-  #           Calculation Methods, in Evaporation and Evapotranspiration: 
+
+  # Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure
+  #           Calculation Methods, in Evaporation and Evapotranspiration:
   #           Measurements and Estimations, Springer, London.
   #             vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea
   #             where:
   #                 tc = average daily air temperature, deg C
   #                 ea  = actual vapor pressure, Pa
-  
+
   kR  = 8.3143   # universal gas constant, J/mol/K (Allen, 1973)
   kMv = 18.02    # molecular weight of water vapor, g/mol (Tsilingiris, 2008)
   kMa = 28.963   # molecular weight of dry air, g/mol (Tsilingiris, 2008)
-  
+
   # calculate the mass mixing ratio of water vapor to dry air (dimensionless)
   wair <- qair / (1 - qair)
-  
-  # calculate water vapor pressure 
+
+  # calculate water vapor pressure
   rv <- kR / kMv
   rd <- kR / kMa
   eact = patm * wair * rv / (rd + wair * rv)
-  
+
   return( eact )
-}
+}

R/calc_vpd.R

@@ -6,15 +6,15 @@
 #' @param qair Air specific humidity (g g-1)
 #' @param eact Water vapour pressure (Pa)
 #' @param tc temperature, deg C
-#' @param tmin (optional) min daily air temp, deg C 
+#' @param tmin (optional) min daily air temp, deg C
 #' @param tmax (optional) max daily air temp, deg C
 #' @param patm Atmospehric pressure (Pa)
-#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing 
+#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing
 #' for calculating it based on standard sea level pressure)
 #'
 #' @return vapor pressure deficit (Pa)
 #' @export
-#' 
+#'
 calc_vpd <- function(
   qair = NA,
   eact = NA,
@@ -24,10 +24,10 @@ calc_vpd <- function(
   patm = NA,
   elv = NA
 ) {
-  
+
   ##-----------------------------------------------------------------------
-  ## Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure 
-  ##           Calculation Methods, in Evaporation and Evapotranspiration: 
+  ## Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure
+  ##           Calculation Methods, in Evaporation and Evapotranspiration:
   ##           Measurements and Estimations, Springer, London.
   ##             vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea
   ##             where:
@@ -37,18 +37,18 @@ calc_vpd <- function(
 
   ## calculate atmopheric pressure (Pa) assuming standard conditions at sea level (elv=0)
   # if (is.na(elv) && is.na(patm) && is.na(eact)){
-  #   
+  #
   #   warning("calc_vpd(): Either patm or elv must be provided if eact is not given.")
   #   vpd <- NA
-  #   
+  #
   # } else {
-    
+
     # if (is.na(eact)){
       patm <- ifelse(is.na(patm),
-                     ingestr::calc_patm(elv),
+                     calc_patm(elv),
                      patm)
     # }
-    
+
     ## Calculate VPD as mean of VPD based on Tmin and VPD based on Tmax if they are availble.
     ## Otherwise, use just tc for calculating VPD.
     vpd <- ifelse(
@@ -62,31 +62,31 @@ calc_vpd <- function(
 }
 
 #' Calculate instantenous VPD from ambient conditions
-#' 
+#'
 #' Follows Abtew and Meleese (2013), Ch. 5 Vapor Pressure Calculation Methods,
 #' in Evaporation and Evapotranspiration
-#' 
+#'
 #' @param qair Air specific humidity (g g-1)
 #' @param eact Water vapour pressure (Pa)
 #' @param tc temperature, deg C
 #' @param patm Atmospehric pressure (Pa)
-#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing 
+#' @param elv Elevation above sea level (m) (Used only if \code{patm} is missing
 #' for calculating it based on standard sea level pressure)
 #'
 #' @return instantenous VPD from ambient conditions
 #' @export
 
 calc_vpd_inst <- function(
   qair=NA,
-  eact=NA, 
+  eact=NA,
   tc=NA,
   patm=NA,
   elv=NA
 ) {
-  
+
   ##-----------------------------------------------------------------------
-  ## Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure 
-  ##           Calculation Methods, in Evaporation and Evapotranspiration: 
+  ## Ref:      Eq. 5.1, Abtew and Meleese (2013), Ch. 5 Vapor Pressure
+  ##           Calculation Methods, in Evaporation and Evapotranspiration:
   ##           Measurements and Estimations, Springer, London.
   ##             vpd = 0.611*exp[ (17.27 tc)/(tc + 237.3) ] - ea
   ##             where:
@@ -101,38 +101,38 @@ calc_vpd_inst <- function(
     calc_eact(qair, patm),
     eact
   )
-  
+
   ## calculate saturation water vapour pressure in Pa
   esat <- 611.0 * exp( (17.27 * tc)/(tc + 237.3) )
-  
+
   ## calculate VPD in units of Pa
-  vpd <- ( esat - eact )    
-  
+  vpd <- ( esat - eact )
+
   ## this empirical equation may lead to negative values for VPD
   ## (happens very rarely). assume positive...
   vpd <- ifelse(
     vpd < 0,
     0,
     vpd
   )
-  
+
   return( vpd )
 }
 
 calc_eact <- function(qair, patm){
-  
+
   # kTo <- 288.15   # base temperature, K (Prentice, unpublished)
   kR  <- 8.3143   # universal gas constant, J/mol/K (Allen, 1973)
   kMv <- 18.02    # molecular weight of water vapor, g/mol (Tsilingiris, 2008)
   kMa <- 28.963   # molecular weight of dry air, g/mol (Tsilingiris, 2008)
-  
+
   ## calculate the mass mixing ratio of water vapor to dry air (dimensionless)
   wair <- qair / (1 - qair)
-  
-  ## calculate water vapor pressure 
+
+  ## calculate water vapor pressure
   rv <- kR / kMv
   rd <- kR / kMa
   eact <- patm * wair * rv / (rd + wair * rv)
-  
+
   return(eact)
-}
+}