R/DRYWETAIR.R

In TrenchR: Tools for Microclimate and Biophysical Ecology

#' @title Properties of Dry Air
#'
#' @description The function calculates several properties of dry air and related characteristics shown as output variables below. The program is based on equations from \insertCite{List1971;textual}{TrenchR} and code implementation from NicheMapR \insertCite{Kearney2017,Kearney2020}{TrenchR}.
#'   \cr \cr
#'  The user must supply values for the input variables \code{db}, \code{bp}, and \code{alt}. If \code{alt} is known (-1000 < alt < 20000) but not BP, then set \code{bp = 0}.
#' 
#' @param db \code{numeric} dry bulb temperature (C).
#' 
#' @param bp \code{numeric} barometric pressure (Pa).
#' 
#' @param alt \code{numeric} altitude (m).
#' 
#' @return Named \code{list} with elements:
#'   \describe{
#'     \item{\code{patmos}}{\code{numeric} standard atmospheric pressure (Pa)}
#'     \item{\code{density}}{\code{numeric} density (\ifelse{html}{\out{kg m<sup>-3</sup>}}{\eqn{kg m^-3}{ASCII}})}
#'     \item{\code{visdyn}}{\code{numeric} dynamic viscosity (\ifelse{html}{\out{kg m<sup>-1</sup> s<sup>-1</sup>}}{\eqn{kg m^-1 s^-1}{ASCII}})}
#'     \item{\code{viskin}}{\code{numeric} kinematic viscosity (\ifelse{html}{\out{m<sup>2</sup> s<sup>-1</sup>}}{\eqn{m^2 s^-1}{ASCII}})}
#'     \item{\code{difvpr}}{\code{numeric} diffusivity of water vapor in air (\ifelse{html}{\out{m<sup>2</sup> s<sup>-1</sup>}}{\eqn{m^2 s^-1}{ASCII}})}
#'     \item{\code{thcond}}{\code{numeric} thermal conductivity (\ifelse{html}{\out{W K<sup>-1</sup> m<sup>-1</sup>}}{\eqn{W K^-1 m^-1}{ASCII}})}
#'     \item{\code{htovpr}}{\code{numeric} latent heat of vaporization of water (\ifelse{html}{\out{J kg<sup>-1</sup>}}{\eqn{J kg^-1}{ASCII}})}
#'     \item{\code{tcoeff}}{\code{numeric} temperature coefficient of volume expansion (\ifelse{html}{\out{K<sup>-1</sup>}}{\eqn{K^-1}{ASCII}})}
#'     \item{\code{ggroup}}{\code{numeric} group of variables in Grashof number (\ifelse{html}{\out{m<sup>-3</sup> K<sup>-1</sup>}}{\eqn{m^-3 K^-1}{ASCII}})}
#'     \item{\code{bbemit}}{\code{numeric} black body emittance (\ifelse{html}{\out{W m<sup>-2</sup>}}{\eqn{W m^-2}{ASCII}})}
#'     \item{\code{emtmax}}{\code{numeric} wave length of maximum emittance (m)}
#'   }
#'
#' @export
#' 
#' @references
#'   \insertAllCited{}
#' 
#' @examples
#'   DRYAIR(db = 30, 
#'          bp = 100*1000, 
#'          alt = 0)
#' 
DRYAIR <- function(db, 
                   bp = 0, 
                   alt = 0) {
  
  stopifnot(alt >= 0, 
            bp  >=  0)

  tstd <- 273.15
  
  pstd   <- 101325.
  patmos <- pstd*((1 - (0.0065 * alt / 288))^(1 / 0.190284))
  
  bp          <- rep(bp, length(patmos))
  bp[bp <= 0] <- patmos[bp <= 0]
  
  density <- bp / (287.04 * (db + tstd))
  
  visnot <- 1.8325E-5
  tnot   <- 296.16
  c      <- 120.
  
  visdyn <- (visnot * ((tnot + c) / (db + tstd + c))) * (((db + tstd) / tnot)^1.5)
  viskin <- visdyn / density
  
  difvpr <- 2.26E-5 * (((db + tstd) / tstd)^1.81) * (1.E5 / bp)
  thcond <- 0.02425 + (7.038E-5 * db)
  htovpr <- 2.5012E6 - 2.3787E3 * db
  tcoeff <- 1. / (db + tstd)
  ggroup <- 0.0980616 * tcoeff / (viskin * viskin)
  bbemit <- 5.670367E-8 * ((db + tstd)^4)
  emtmax <- 2.897E-3 / (db + tstd)
  
  list(patmos  = patmos, 
       density = density, 
       visdyn  = visdyn, 
       viskin  = viskin, 
       difvpr  = difvpr, 
       thcond  = thcond, 
       htovpr  = htovpr, 
       tcoeff  = tcoeff, 
       ggroup  = ggroup, 
       bbemit  = bbemit, 
       emtmax  = emtmax)

}


#' @title Saturation Vapor Pressure
#'
#' @description The function calculates saturation vapor pressure for a given air temperature. The program is based on equations from \insertCite{List1971;textual}{TrenchR} and code implementation from NicheMapR \insertCite{Kearney2017,Kearney2020}{TrenchR}.
#' 
#' @param db \code{numeric} Dry bulb temperature (C)
#' 
#' @return esat \code{numeric} Saturation vapor pressure (Pa)
#' 
#' @export
#' 
#' @references
#'   \insertAllCited{}
#' 
#' @examples
#'   VAPPRS(db = 30)
#' 
VAPPRS <- function(db) {
  
  t <- celsius_to_kelvin(db)
  
  loge <- t
  loge[t <= 273.16] <- -9.09718 * (273.16 / t[t <= 273.16] - 1.) - 3.56654 * log10(273.16 / t[t <= 273.16]) + 0.876793 * (1. - t[t <= 273.16] / 273.16) + log10(6.1071)
  loge[t > 273.16]  <- -7.90298 * (373.16 / t[t > 273.16] - 1.) + 5.02808 * log10(373.16 / t[t > 273.16]) - 1.3816E-07 * (10.^(11.344 * (1. - t[t > 273.16] / 373.16)) - 1.) + 8.1328E-03 * (10.^(-3.49149 * (373.16 / t[t > 273.16] - 1.)) - 1.) + log10(1013.246)
  
  (10.^loge) * 100
  
}

#' @title Properties of Wet Air
#'
#' @description The function calculates several properties of humid air described as output variables below. The program is based on equations from \insertCite{List1971;textual}{TrenchR} and code implementation from NicheMapR \insertCite{Kearney2017,Kearney2020}{TrenchR}.
#'  \cr 
#'  WETAIR must be used in conjunction with \code{\link{VAPPRS}}. Input variables are shown below. See Details. 
#'
#' @details The user must supply known values for DB and BP (BP at one standard
#'   atmosphere is 101,325 pascals). Values for the remaining variables are
#'   determined by whether the user has either (1) psychrometric data (WB or RH),
#'   or (2) hygrometric data (DP):
#'   \describe{
#'     \item{Psychrometric data}{If WB is known but not RH, then set RH = -1
#'       and DP = 999. If RH is known but not WB then set WB = 0 and DP = 999.}
#'     \item{Hygrometric data}{If DP is known, set WB = 0 and RH = 0.}
#'   }
#'
#' @param db \code{numeric} dry bulb temperature (C).
#' 
#' @param wb \code{numeric} wet bulb temperature (C).
#' 
#' @param rh \code{numeric} relative humidity (\%).
#' 
#' @param dp \code{numeric} dew point temperature (C).
#' 
#' @param bp \code{numeric} barometric pressure (Pa).
#' 
#' @return Named \code{list} with elements:
#'   \describe{
#'     \item{\code{e}}{\code{numeric} saturation vapor pressure (Pa)}
#'     \item{\code{vd}}{\code{numeric} vapor density (\ifelse{html}{\out{kg m<sup>-3</sup>}}{\eqn{kg m^-3}{ASCII}})}
#'     \item{\code{rw}}{\code{numeric} mixing ratio (\ifelse{html}{\out{kg kg<sup>-1</sup>}}{\eqn{kg kg^-1}{ASCII}})}
#'     \item{\code{tvir}}{\code{numeric} virtual temperature (K)}
#'     \item{\code{tvinc}}{\code{numeric} virtual temperature increment (K)}
#'     \item{\code{denair}}{\code{numeric} density of the air (\ifelse{html}{\out{kg m<sup>-3</sup>}}{\eqn{kg m^-3}{ASCII}})}
#'     \item{\code{cp}}{\code{numeric} specific heat of air at constant pressure (\ifelse{html}{\out{J kg<sup>-1</sup> K<sup>-1</sup>}}{\eqn{J kg^-1 K^-1}{ASCII}})}
#'     \item{\code{wtrpot}}{\code{numeric} water potential (Pa)}
#'     \item{\code{rh}}{\code{numeric} relative humidity (\%)}
#'   }
#' 
#' @export
#' 
#' @references
#'   \insertAllCited{}
#' 
#' @examples
#'   WETAIR(db = 30, 
#'          wb = 28, 
#'          rh = 60, 
#'          bp = 100 * 1000)
#' 
WETAIR <- function (db, 
                    wb = db, 
                    rh = 0, 
                    dp = 999, 
                    bp = 101325) {
  
  stopifnot(rh <= 100, 
            bp >  0)
  
  tk    <- celsius_to_kelvin(db)
  esat  <- VAPPRS(db)

  if(dp < 999.0) {
    
    e  <- VAPPRS(dp)
    rh <- (e / esat) * 100
    
  } else {
    
    if(min(rh) > -1){
      
      e <- esat * rh / 100
      
    } else{
      
      wbd   <- db - wb
      wbsat <- VAPPRS(wb)
      dltae <- 0.000660 * (1.0 + 0.00115 * wb) * bp * wbd
      e     <- wbsat - dltae
      rh    <- (e / esat) * 100
      
    }
  }
  
  rw     <- ((0.62197 * 1.0053 * e) / (bp - 1.0053 * e))
  vd     <- e * 0.018016 / (0.998 * 8.31434 * tk)
  tvir   <- tk * ((1.0 + rw / (18.016 / 28.966)) / (1.0 + rw))
  tvinc  <- tvir - tk
  denair <- 0.0034838 * bp / (0.999 * tvir)
  cp     <- (1004.84 + (rw * 1846.40)) / (1.0 + rw)
  
  if (min(rh) <= 0.0) {
    
    wtrpot <- -999
    
  } else {
    
    wtrpot <- 4.615e+5 * tk * log(rh / 100.0)
    
  }
  
  list(e      = e, 
       esat   = esat, 
       vd     = vd, 
       rw     = rw, 
       tvinc  = tvinc, 
       denair = denair, 
       cp     = cp, 
       wtrpot = wtrpot, 
       rh     = rh)
  
}