From b503a82329a6069033d49b7eb7ee69e4aa26bea3 Mon Sep 17 00:00:00 2001
From: Dongdong Kong <kongdd@users.noreply.github.com>
Date: Mon, 25 Sep 2023 14:08:59 +0800
Subject: [PATCH] fix the error rho_a

---
 R/PML.R                  | 10 +++++-----
 R/constant.R             |  3 ++-
 R/tidy_forcing_flux.R    |  6 +++---
 man/PML.Rd               |  2 +-
 man/PML_calib.Rd         |  2 +-
 man/tidy_forcing_flux.Rd |  2 +-
 6 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/R/PML.R b/R/PML.R
index ce65c34..b8503e6 100644
--- a/R/PML.R
+++ b/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]
diff --git a/R/constant.R b/R/constant.R
index 65c8535..708064f 100644
--- a/R/constant.R
+++ b/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
diff --git a/R/tidy_forcing_flux.R b/R/tidy_forcing_flux.R
index 6d0225d..f1d412b 100644
--- a/R/tidy_forcing_flux.R
+++ b/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)
 }
diff --git a/man/PML.Rd b/man/PML.Rd
index 67283c7..3844ff0 100644
--- a/man/PML.Rd
+++ b/man/PML.Rd
@@ -52,7 +52,7 @@ parameter.
 (optional, for speeding up optimization)
 \itemize{
 \item \code{lambda}  :
-\item \code{rou_a}   :
+\item \code{rho_a}   : air density, \verb{[g m-3]}
 \item \code{gama}    :
 \item \code{epsilon} :
 }}
diff --git a/man/PML_calib.Rd b/man/PML_calib.Rd
index de25a7c..723045e 100644
--- a/man/PML_calib.Rd
+++ b/man/PML_calib.Rd
@@ -50,7 +50,7 @@ parameter.
 (optional, for speeding up optimization)
 \itemize{
 \item \code{lambda}  :
-\item \code{rou_a}   :
+\item \code{rho_a}   : air density, \verb{[g m-3]}
 \item \code{gama}    :
 \item \code{epsilon} :
 }}
diff --git a/man/tidy_forcing_flux.Rd b/man/tidy_forcing_flux.Rd
index c085dfc..96b1e55 100644
--- a/man/tidy_forcing_flux.Rd
+++ b/man/tidy_forcing_flux.Rd
@@ -24,7 +24,7 @@ tidy_forcing_flux(df)
 \item lambda
 \item Rn
 \item Eeq
-\item rou_a
+\item rho_a
 \item gama
 \item epsilon
 \item fval_soil