diff --git a/functions/.Rapp.history b/functions/.Rapp.history deleted file mode 100644 index 895abff..0000000 --- a/functions/.Rapp.history +++ /dev/null @@ -1,166 +0,0 @@ -calc_km_pa = function(temp, z) {# - patm = calc_patm(z) # - rat = patm / calc_patm(0)# - R = 8.314 # - O2 = 2.09476e5 # - Kc25 = 41.03 * rat # - Ko25 = 28210 * rat # - Hkc = 79430 # - Hko = 36380 # - temp_k = 273.15 + temp# -# - Kc_pa =Kc25 * exp(Hkc * ((temp_k - 298.15) / (298.15 * R * temp_k)))# - Ko_pa =Ko25* exp(Hko * ((temp_k - 298.15) / (298.15 * R * temp_k)))# - O2_pa = O2 * (1e-6) * patm # - Km_pa = Kc_pa * (1 + O2_pa/Ko_pa)# - Km_pa # -} -calc_km_pa(25, 0) -calc_patm = function(z) {# - kPo = 101325 # standard atmosphere, Pa (Allen, 1973)# - kTo = 298.15 # base temperature, K (Prentice, unpublished)# - kL = 0.0065 # temperature lapse rate, K/m (Allen, 1973)# - kG = 9.80665 # gravitational acceleration, m/s**2 (Allen, 1973)# - kR = 8.3143 # universal gas constant, J/mol/K (Allen, 1973)# - kMa = 0.028963 # molecular weight of dry air, kg/mol (Tsilingiris, 2008)# - patm = kPo*(1.0 - kL*z/kTo)**(kG*kMa/(kR*kL))# - patm# -} -calc_km_pa(25, 0) -log(146) + log(71.9) - log (1.6) -(log(146) + log(71.9) - log (1.6))/ 2 -(log(146) + log(.0719) - log (1.6))/ 2 -(log(240) + log(.0719) - log (1.6))/ 2 -calc_chi = function(temp, z, vpdo, cao){ # temp in °C, z in m, vpd in kPa, ca in ppm# - beta = 30 # - patm = calc_patm(z) # - vpd = calc_vpd(temp, z, vpdo) # - vpd_pa = vpd * 1000# - ca_pa = cao * 1e-6 * patm# - km = calc_km_pa(temp, z) # - gammastar = calc_gammastar_pa(temp, z)# - nstar = calc_nstar(temp, z)# - xi = sqrt((beta * (km + gammastar)) / (1.6 * nstar))# - chi = (gammastar / ca_pa) + (1 - (gammastar / ca_pa)) * (xi / (xi + sqrt(vpd_pa)))# - chi # -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_vpd = function(temp, z, vpdo){# - vpdo_pa = vpdo * 1000 # VPD at sea level (assumed to be the case for CRU)# - po = 101325 # Atmospheric pressure under standard conditions (i.e., sea level)# - patm = calc_patm(z) # actual atmospheric pressure# - es = 610.8 * exp((17.27 * temp) / (237.3 + temp)) # saturation vapor pressure (Pa) http://cronklab.wikidot.com/calculation-of-vapour-pressure-deficit# - eao = -(vpdo_pa - es) # actual vapor pressure at sea level (Pa)# - vpd = es - eao * (patm / po) # VPD at z (Pa)# - vpd / 1000 # VPD at z (kPa)# -# -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_gammastar_pa = function(temp, z) {# - patm = calc_patm(z)# - rat = calc_patm(z) / calc_patm(0)# - #gammastar25 = 42.75 # ppm# - gammastar25 = 4.332 * rat # Pa# - Hgm=37830 # J mol-1# - R = 8.314 # J K-1 mol-1# - O2 = 2.09476e5 # ppm# - O2_0 = O2 * 1e-6 * calc_patm(0)# - O2_z = O2 * 1e-6 * calc_patm(z)# - temp_k = 273.15+ temp# - gStar_pa = gammastar25*exp((Hgm/R)*(1/298.15-1/temp_k))# - #gStar_ppm * 1e-6 * patm #convert to Pa# - gStar_pa * (O2_z / O2_0) # vary based on oxygen due to Rubisco specificity factor# -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_nstar = function(temp, z){ # temp in °C and z in m# - patm = calc_patm(z)# - # viscosity correction factor = viscosity( temp, press )/viscosity( 25 degC, 1013.25 Pa) # - ns = calc_viscosity_h2o( temp, z ) # Pa s # - ns25 = calc_viscosity_h2o( 25, 1013.25 ) # Pa s # - nstar = ns / ns25 # (unitless)# - nstar# -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_viscosity_h2o = function(temp, z){ #temp in °C and z in m# - tk_ast = 647.096 # Kelvin# - rho_ast = 322.0 # kg/m**3# - mu_ast = 1e-6 # Pa s# - patm = calc_patm(z)# - rho = calc_density_h2o(temp, z) # density of water (kg m-3)# - # Calculate dimensionless parameters:# - tbar = (temp + 273.15)/tk_ast# - tbarx = tbar**(0.5)# - tbar2 = tbar**2# - tbar3 = tbar**3# - rbar = rho/rho_ast# - # Calculate mu0 (Eq. 11 & Table 2, Huber et al., 2009):# - mu0 = 1.67752 + 2.20462/tbar + 0.6366564/tbar2 - 0.241605/tbar3# - mu0 = 1e2*tbarx/mu0# - # Create Table 3, Huber et al. (2009):# -# - h_array1 = c(0.520094, 0.0850895, -1.08374, -0.289555, 0.0, 0.0)# - h_array2 = c(0.222531, 0.999115, 1.88797, 1.26613, 0.0, 0.120573)# - h_array3 = c(-0.281378, -0.906851, -0.772479, -0.489837, -0.257040, 0.0)# - h_array4 = c(0.161913, 0.257399, 0.0, 0.0, 0.0, 0.0)# - h_array5 = c(-0.0325372, 0.0, 0.0, 0.0698452, 0.0, 0.0)# - h_array6 = c(0.0, 0.0, 0.0, 0.0, 0.00872102, 0.0)# - h_array7 = c(0.0, 0.0, 0.0, -0.00435673, 0.0, -0.000593264)# - h_array = rbind(h_array1, h_array2, h_array3, h_array4, h_array5, h_array6, h_array7)# - # Calculate mu1 (Eq. 12 & Table 3, Huber et al., 2009):# - mu1 = 0.0# - ctbar = (1.0/tbar) - 1.0# - for (i in 1:6){# - coef1 = ctbar**(i-1)# - coef2 = 0.0# - for (j in 1:7){# - coef2 = coef2 + h_array[j,i] * (rbar - 1.0)**(j-1)# - }# - mu1 = mu1 + coef1 * coef2 # - }# - mu1 = exp( rbar * mu1 )# -# - # Calculate mu_bar (Eq. 2, Huber et al., 2009)# - # assumes mu2 = 1# - mu_bar = mu0 * mu1# -# - # Calculate mu (Eq. 1, Huber et al., 2009)# - viscosity_h2o = mu_bar * mu_ast # Pa s# - viscosity_h2o# -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_density_h2o = function(temp, z){ # temp in °C and z in m# - patm = calc_patm(z)# - # Calculate lambda, (bar cm**3)/g:# - my_lambda = 1788.316 + 1.55053*temp + -0.4695911*temp*temp + (3.096363e-3)*temp*temp*temp + -(7.341182e-6)*temp*temp*temp*temp# -# - # Calculate po, bar# - po = 5918.499 + 58.05267*temp + -1.1253317*temp*temp + (6.6123869e-3)*temp*temp*temp + -(1.4661625e-5)*temp*temp*temp*temp# -# - # Calculate vinf, cm**3/g# - vinf = 0.6980547 + -(7.435626e-4)*temp + (3.704258e-5)*temp*temp + -(6.315724e-7)*temp*temp*temp + (9.829576e-9)*temp*temp*temp*temp + -(1.197269e-10)*temp*temp*temp*temp*temp + (1.005461e-12)*temp*temp*temp*temp*temp*temp + -(5.437898e-15)*temp*temp*temp*temp*temp*temp*temp + (1.69946e-17)*temp*temp*temp*temp*temp*temp*temp*temp + -(2.295063e-20)*temp*temp*temp*temp*temp*temp*temp*temp*temp# -# - # Convert pressure to bars (1 bar = 100000 Pa)# - pbar = (1e-5)*patm# - # Calculate the specific volume (cm**3 g**-1):# - vau = vinf + my_lambda/(po + pbar)# -# - # Convert to density (g cm**-3) -> 1000 g/kg; 1000000 cm**3/m**3 -> kg/m**3:# - density_h2o = (1e3/vau)# - density_h2o# -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -calc_chi = function(temp, z, vpdo, cao){ # temp in °C, z in m, vpd in kPa, ca in ppm# - beta = 146 # - patm = calc_patm(z) # - vpd = calc_vpd(temp, z, vpdo) # - vpd_pa = vpd * 1000# - ca_pa = cao * 1e-6 * patm# - km = calc_km_pa(temp, z) # - gammastar = calc_gammastar_pa(temp, z)# - nstar = calc_nstar(temp, z)# - xi = sqrt((beta * (km + gammastar)) / (1.6 * nstar))# - chi = (gammastar / ca_pa) + (1 - (gammastar / ca_pa)) * (xi / (xi + sqrt(vpd_pa)))# - chi # -} -calc_chi(temp = 25, z = 0, vpdo =1, cao=400) -source('calc_optimal_vcmax.R')