fix the error rho_a

R/PML.R

@@ -27,7 +27,7 @@
 #'
 #' (optional, for speeding up optimization)
 #' - `lambda`  :
-#' - `rou_a`   :
+#' - `rho_a`   : air density, `[g m-3]`
 #' - `gama`    :
 #' - `epsilon` :
 #'
@@ -146,13 +146,13 @@ PML <- function(
   lambda <- data$lambda
   if (is.null(lambda)) lambda <- get_lambda(Tavg) # [2500 J g-1]
 
-  rou_a <- data$rou_a
+  rho_a <- data$rho_a
   gama <- data$gama
   epsilon <- data$epsilon
   Eeq <- data$Eeq
 
-  if (is.null(rou_a)) {
-    rou_a   <- 3.846 * 10^3 * Pa / (Tavg + 273.15)
+  if (is.null(rho_a)) {
+    rho_a   <- 3.486 * 10^3 * Pa / (1.01 * (Tair + 273)) # g m-3
     gama    <- Cp * Pa / (0.622 * lambda)
     slop    <- cal_slope(Tavg) # kPa deg-1
     epsilon <- slop / gama
@@ -227,7 +227,7 @@ PML <- function(
   # Transpiration from plant cause by radiation water transfer
   LEcr <- epsilon * Rn * (1 - Tou) / (epsilon + 1 + Ga / Gc) # W m-2
   # Transpiration from plant cause by aerodynamic water transfer
-  LEca <- (rou_a * Cp * Ga * VPD / gama) / (epsilon + 1 + Ga / Gc) # W m-2
+  LEca <- (rho_a * Cp * Ga * VPD / gama) / (epsilon + 1 + Ga / Gc) # W m-2
 
   Ecr  <- LEcr / lambda * 86400 * 10^-3 # [W m-2] change to [mm d-1]
   Eca  <- LEca / lambda * 86400 * 10^-3 # [W m-2] change to [mm d-1]

R/constant.R

@@ -2,4 +2,5 @@
 # derivied from Kelliher FM et al. (1995), AFM
 kmar <- 0.40        # von Karman's constant 0.40
 Zob  <- 15          # m, making sure higher than hc
-Cp   <- 4.2 * 0.242 # specific heat at constant pressure, 1.013/1000  [MJ/kg/degC]
+# Cp   <- 4.2 * 0.242 # specific heat at constant pressure, 1.013/1000  [MJ/kg/degC]
+Cp   <- 1.013       # [J g-1 degC-1], v20230925

R/tidy_forcing_flux.R

@@ -15,7 +15,7 @@
 #' - lambda
 #' - Rn
 #' - Eeq
-#' - rou_a
+#' - rho_a
 #' - gama
 #' - epsilon
 #' - fval_soil
@@ -47,7 +47,7 @@ tidy_forcing_flux <- function(df) # , par
   Rn <- get_Rn(df$Rs, df$Rl_in, Tavg, df$Albedo, df$Emiss)
 
   # 2. intermediate variables
-  rou_a <- 3.846 * 10^3 * Pa / (Tavg + 273.15) # kg m-3
+  rho_a <- 3.486 * 10^3 * Pa / (1.01 * (Tair + 273)) # g m-3
   gama <- Cp * Pa / (0.622 * lambda) # kpa deg-1
   slop <- 4098 * 0.6108 * exp((17.27 * Tavg) / (Tavg + 237.3)) / (Tavg + 237.3)^2 # kpa deg-1
   epsilon <- slop / gama
@@ -69,6 +69,6 @@ tidy_forcing_flux <- function(df) # , par
   #     G_soil       = (Tavg-Tavg_1) * 1/0.408*0.38 / (lambda * 1e-3);# MJ /m2/day
   #     G_soil       = G_soil / (1 / lambda*86400*10^-3);# W/m2
   # }
-  out <- data.table(ETobs, GPPobs, Rn, Eeq, lambda, rou_a, gama, epsilon, fval_soil)
+  out <- data.table(ETobs, GPPobs, Rn, Eeq, lambda, rho_a, gama, epsilon, fval_soil)
   cbind(df, out)
 }