R/get.PAR.R

In gholtgrieve/gassyPants: Functions Useful for Working With Gasses Dissolved in Water

#' @title Estimates Photosynthetically Active Radiation at the Earth's Surface
#' @description
#' Estimates 'blue-sky' photosynthetically active radiation (PAR) at the Earth's surface
#' for a given location and timeframe. Considers both direct and diffuse radiation,
#' but does not consider topographic or atmospheric effects.\cr
#' @usage
#' get.PAR(lat, long, localHour, DOY, tz, masl = 0.001, aspect = 0.001, slope = 0.001, trans = 0.7)
#' @param lat Latitude in decimal degrees
#' @param long Longitude in decimal degrees
#' @param localHour Local time of day expressed as a decimal (i.e., 1:15 pm == 13.25)
#' @param DOY Ordinal day of year ranging between 1 and 366 (1 == Jan 1; 366 == Dec 31 on a leap year.)
#' @param tz Time zone expressed as hours from UTC (i.e., Pacific Standard Time == -8)
#' @param masl Meters above sea level
#' @param aspect Aspect of the land being modeled in degrees. Defaults to 180 (south).
#' @param slope Slope of the land being modeled in degree. Defaults to 0.001 (flat).
#' @param trans Atmospheric transmissivity (Fraction of extra-terrestrial solar radiation reaching the Earth's surface under cloudless
#' conditions). Ranges from 0-1.  Defaults to 0.7, which is a very approximate average.
#' @return Photosynthetically active radiation (PAR) in units of microEinstiens s-1 m-2
#' @note Specific PAR values are aproximations based on ideal (and probably unrealistic) condtions.
#' @author Gordon W. Holtgrieve (gholt@uw.edu)
#' @references
#' Venkiteswaran et al. 2007
#' Gates 1966
#' Jellison and Melack 1993
#' McCree 1972
#' @export

#*****PAR*********************************************
#Estimates the amount of light at surface of water body. Units of microEinstiens s-1 m-2.
#**********************************************************************
# Created by GWH on 8/13/2010
# Based on BaMM light model, which is based on Venkiteswaran et al. 2007

# SolConst              'Solar constant in Watts m-2
# SolarHour             'solar hour for a given position on the globe
# HourAngle As Double   'movement of sun around earth in degrees
# SolDec As Double      'annual fluctuation of sun between two poles, varies between -23.45 and +23.45 radians
# SolAlt As Double      'height of sun above horizon in radians
# SolAzi As Double      'angle in degrees between solar beams and North-South axis of Earth
# SolAziCorr As Double  'corrected solar azimuth
# CosInci As Double     'angle between incoming solar beams and aspect and slope of the surface
# OptCorr As Double     'correct the relative path length of the optical air mass for masl
# AtmPressCorr As Double 'atmospheric pressure at a given masl
# SOuter As Double      'no idea
# SNormal As Double     'radiation normal to the beam
# SDirect As Double     'direct radation on earth's surface corrected for orientation and slope of surface
# SDiffuse As Double    'Indirect or diffuse sunlight for each solar altitude
# STotalWatts As Double 'total radiation to the water body in watts m-2
# PAR As Double         'total radiation to the water body in microEinsteins s-1 m-2

get.PAR <- function (lat, long, localHour, DOY, tz, masl=0.001,
                     aspect=0.001, slope=0.001, trans=0.7) {

  # SolConst      solar constant in Watts m-2
  # SolarHour     solar hour for a given position on the globe
  # HourAngle     movement of sun around earth in degrees
  # SolDec        annual fluctuation of sun between two poles, varies between -23.45 and +23.45 radians
  # SolAlt        height of sun above horizon in radians
  # SolAzi        angle in degrees between solar beams and North-South axis of Earth
  # SolAziCorr    corrected solar azimuth
  # CosInci       angle between incoming solar beams and aspect and slope of the surface
  # OptCorr       correct the relative path length of the optical air mass for masl
  # AtmPressCorr  atmospheric pressure at a given masl
  # SOuter        no idea
  # SNormal       radiation normal to the beam
  # SDirect       direct radation on earth's surface corrected for orientation and slope of surface
  # SDiffuse      Indirect or diffuse sunlight for each solar altitude
  # STotalWatts   total radiation to the water body in watts m-2
  # PAR           total radiation to the water body in microEinsteins s-1 m-2

  SolConst <- 1367

  SolarHour <- localHour + 4 / 60 * long - tz
  HourAngle <- 15 * (SolarHour - 12)
  SolDec <- -23.45 * cos(360 * (DOY + 10) / 365 * pi/180)
  SolAlt <- asin(sin(lat * pi/180) * sin(SolDec * pi/180) + cos(lat * pi/180) * cos(SolDec * pi/180) * cos(HourAngle * pi/180))
  SolAzi <- (sin(SolDec * pi/180) * cos(lat * pi/180) * sin(lat * pi/180) * cos(HourAngle * pi/180)) / cos(SolAlt)

  SolAzi[SolAzi > 0.999] <- 0.999
  SolAzi <- SolAzi / pi/ 180

  SolAziCorr <- 360 - SolAzi
  SolAziCorr[HourAngle < 12] <- SolAzi[HourAngle < 12]

  CosInci <- sin(SolAlt) * cos(slope * pi/180) + cos(SolAlt) * sin(slope * pi/180) * cos((SolAziCorr - aspect) * pi/180)
  AtmPressCorr <- ((288.15 - 0.0065 * masl) / 288.15) ^ 5.256
  OptCorr <- AtmPressCorr * trans ^ ((1229 + (614 * sin(SolAlt)) ^ 2) ^ 0.5 - 614 * sin(SolAlt))
  SOuter <- SolConst * (1 + 0.034 * cos(2 * pi * DOY / 365))
  SNormal <- SOuter * OptCorr

  SDirect <- SNormal * CosInci
  SDirect[SDirect < 0] <- 0

  SDiffuse <- SOuter * (0.271 - 0.294 * OptCorr) * sin(SolAlt)
  SDiffuse[SDiffuse < 0] <- 0

  STotalWatts <- SDirect + SDiffuse     #Full spectrun irradiance (FSI) in watts m-2
# Convert FSI to photosynthetically active radiation
# (PAR=0.45 X FSI; Gates 1966, Jellison and Melack 1993).
# Convert to microEinsteins s-1 m-2.
# Rough conversion is 1 W m-2 (PAR) <-> 4.57 uE m-2 s-1 (PAR) (McCree 1972).
  PAR <- STotalWatts * 0.45 * 4.57
  return(PAR)
}