Addition of home temperature to Biomee and P-model

.gitignore

@@ -6,6 +6,13 @@
 .DS_Store
 .idea
 
+# VS Code local settings
+.vscode/
+
+# act local config & secrets
+.actrc
+.secrets
+
 # source files
 src/*.o
 src/*.so
@@ -18,4 +25,5 @@ vignettes/*.html
 # documentation dir
 docs/
 
+# CRAN submission archives
 CRAN-SUBMISSION

DESCRIPTION

@@ -66,6 +66,7 @@ Imports:
     magrittr,
     GenSA,
     BayesianTools,
+    lubridate,
     multidplyr,
     stats,
     utils

NAMESPACE

@@ -14,6 +14,7 @@ export(runread_pmodel_f)
 import(BayesianTools)
 import(GenSA)
 import(dplyr)
+import(lubridate)
 importFrom(magrittr,"%>%")
 importFrom(stats,setNames)
 useDynLib(rsofun)

NEWS.md

@@ -5,6 +5,10 @@
   easier processing within cost-functions. Default cost-functions were updated,
   but currently ignore the names. This is fully backward compatible, but allows
   to use the names in user-created cost functions.
+* Internally, P-model now computes temperature dependencies of jmax and vcmax
+  with Kumarathunge et al. (2019) instead of Kattge & Knorr (2007), effectively 
+  using a dampened signal for tc_growth and long-term averages for tc_home. This 
+  alters the simulation results. Formattting of inputs remains unchanged.
 
 # rsofun 5.1.0
 

R/run_biomee_f_bysite.R

@@ -16,6 +16,8 @@
 #' @param luc_forcing An array of land use change (LUC) used during transient phase.
 #'
 #' For further specifications of above inputs and examples see \code{\link{biomee_gs_leuning_drivers}}, \code{\link{biomee_p_model_drivers}}, or \code{\link{biomee_p_model_luluc_drivers}}.
+
+#' @import lubridate 
 #' 
 #' @returns A data.frame with columns containing model output for each land unit (LU). 
 #' See examples \code{\link{biomee_gs_leuning_output}}, \code{\link{biomee_p_model_output}}, or \code{\link{biomee_p_model_luluc_output}}. 
@@ -231,6 +233,48 @@ run_biomee_f_bysite <- function(
   makecheck = TRUE
 ){
   ndayyear <- 365
+
+  # Default value for tc_home
+  if ("tc_home" %in% names(site_info)) {
+    stop("Unexpectedly received site_info$tc_home; it should be calculated internally.")
+  }
+
+  conditionally_add_tmax <- function(df){
+    need_to_add_tmax <- !("tmax" %in% colnames(df))
+    if (need_to_add_tmax) {
+      df %>% dplyr::group_by(.data$date) %>%
+        dplyr::summarise(daily_tmax = max(.data$temp, na.rm = TRUE)) %>%
+        dplyr::ungroup()
+    } else {
+      df %>% dplyr::rename(daily_tmax = "tmax")
+    }
+  }
+  tc_home <- forcing %>%
+    # conditionally add daily max temp (if needed, e.g. when running "gs_leuning" with hourly forcing)
+    conditionally_add_tmax() %>%
+    # add grouping variables:
+    mutate(month = lubridate::month(.data$date), year = lubridate::year(.data$date)) %>%
+    # monthly means of daily maximum:
+    group_by(.data$year, .data$month) %>%
+    summarise(monthly_avg_daily_tmax = mean(.data$daily_tmax, na.rm = TRUE), .groups = "drop") %>%
+    # warmest month of each year:
+    group_by(year) %>%
+    summarise(t_warmest_month = max(.data$monthly_avg_daily_tmax)) %>%
+    # mean of yearly warmest months:
+    ungroup() %>%
+    summarise(tc_home = mean(.data$t_warmest_month, na.rm = TRUE)) %>%
+    # extract scalar value
+    dplyr::pull(.data$tc_home)
+
+    site_info$tc_home <- tc_home     # TODO: rather in site_info or in params_tile?
+    # params_tile$tc_home <- tc_home # TODO: rather in site_info or in params_tile?
+
+  # Validate calculation
+  if (is.na(site_info$tc_home) || length(site_info$tc_home) == 0) {
+    warning("Calculated tc_home is NA or missing; defaulting to 25C.")
+    site_info$tc_home <- 25
+  }
+
   # record number of years in forcing data
   # frame to use as default values (unless provided othrwise as params_siml$nyeartrend)
   forcing_years <- nrow(forcing)/(ndayyear * params_siml$steps_per_day)
@@ -548,7 +592,8 @@ prepare_site_info <- function(site_info){
   site_info <- site_info %>% select(
     "lon",
     "lat",
-    "elv"
+    "elv",
+    "tc_home"
   )
   return(site_info)
 }

R/run_pmodel_f_bysite.R

@@ -15,6 +15,7 @@
 #' For further specifications of above inputs and examples see \code{\link{p_model_drivers}} or \code{\link{p_model_drivers_vcmax25}}
 
 #' @import dplyr
+#' @import lubridate
 #' 
 #' @returns Model output is provided as a tidy dataframe, with columns:
 #' \describe{
@@ -159,7 +160,46 @@ run_pmodel_f_bysite <- function(
   secs_per_tstep <- difftime(times[1], times[2], units = "secs") %>%
     as.integer() %>%
     abs()
-
+
+  # Default value for tc_home
+  if ("tc_home" %in% names(site_info)) {
+    stop("Unexpectedly received site_info$tc_home; it should be calculated internally.")
+  }
+
+  conditionally_add_tmax <- function(df){
+    need_to_add_tmax <- !("tmax" %in% colnames(df))
+    if (need_to_add_tmax) {
+      df %>% dplyr::group_by(.data$date) %>%
+        dplyr::summarise(daily_tmax = max(.data$temp, na.rm = TRUE)) %>%
+        dplyr::ungroup()
+    } else {
+      df %>% dplyr::rename(daily_tmax = "tmax")
+    }
+  }
+  tc_home <- forcing %>%
+    # conditionally add daily max temp (if needed, e.g. when running "gs_leuning" with hourly forcing)
+    conditionally_add_tmax() %>%
+    # add grouping variables:
+    mutate(month = lubridate::month(.data$date), year = lubridate::year(.data$date)) %>%
+    # monthly means of daily maximum:
+    group_by(.data$year, .data$month) %>%
+    summarise(monthly_avg_daily_tmax = mean(.data$daily_tmax, na.rm = TRUE), .groups = "drop") %>%
+    # warmest month of each year:
+    group_by(year) %>%
+    summarise(t_warmest_month = max(.data$monthly_avg_daily_tmax)) %>%
+    # mean of yearly warmest months:
+    ungroup() %>%
+    summarise(tc_home = mean(.data$t_warmest_month, na.rm = TRUE)) %>%
+    # extract scalar value
+    dplyr::pull(.data$tc_home)
+
+    site_info$tc_home <- tc_home
+  # Validate calculation
+  if (is.na(site_info$tc_home) || length(site_info$tc_home) == 0) {
+    warning("Calculated tc_home is NA or missing; defaulting to 25C.")
+    site_info$tc_home <- 25
+  }
+
   # re-define units and naming of forcing dataframe
   # keep the order of columns - it's critical for Fortran (reading by column number)
   forcing_features <- c(
@@ -307,6 +347,7 @@ run_pmodel_f_bysite <- function(
       latitude                  = as.numeric(site_info$lat),
       altitude                  = as.numeric(site_info$elv),
       whc                       = as.numeric(site_info$whc),
+      tc_home                   = as.numeric(site_info$tc_home),
       n                         = as.integer(nrow(forcing)), # number of rows in matrix (pre-allocation of memory)
       par                       = c(as.numeric(params_modl$kphio), # model parameters as vector in order
                                     as.numeric(params_modl$kphio_par_a),

src/biosphere_pmodel.mod.f90

@@ -136,7 +136,8 @@ function biosphere_annual() result( out_biosphere )
                   myinterface%climate(doy), &
                   myinterface%grid, &
                   init_daily, &
-                  myinterface%params_siml%in_ppfd &
+                  myinterface%params_siml%in_ppfd, &
+                  myinterface%tc_home &
                   )
         ! if (verbose) print*,'... done'
 

src/gpp_biomee.mod.f90

@@ -68,7 +68,7 @@ subroutine gpp( forcing, vegn )
     integer :: i
     real   :: rad_top                                      ! downward radiation at the top of the canopy, W/m2
     real   :: rad_net                                      ! net radiation absorbed by the canopy, W/m2
-    real   :: TairC, TairK                                 ! air temperature, degC and degK
+    real   :: TairC, TairK                                 ! air temperature, deg C and deg K
     real   :: cana_q                                       ! specific humidity in canopy air space, kg/kg
     real   :: cana_co2                                     ! co2 concentration in canopy air space, mol CO2/mol dry air
     real   :: p_surf                                       ! surface pressure, Pa
@@ -149,7 +149,7 @@ subroutine gpp( forcing, vegn )
           rad_net  = f_light(cc%layer) * forcing%radiation * 0.9 ! net radiation absorbed by the canopy, W/m2
           p_surf   = forcing%P_air  ! Pa
           TairK    = forcing%TairK ! K
-          TairC    = forcing%TairC ! degC
+          TairC    = forcing%TairC ! deg C
           cana_q   = (calc_esat(TairC) * forcing%RH * h2o_molmass) / (p_surf * kMa)  ! air specific humidity, kg/kg
           cana_co2 = forcing%CO2 ! co2 concentration in canopy air space, mol CO2/mol dry air
 
@@ -240,6 +240,7 @@ subroutine gpp( forcing, vegn )
                                 tc             = vegn%dampended_forcing%temp, &
                                 vpd            = vegn%dampended_forcing%vpd, &
                                 patm           = vegn%dampended_forcing%patm, &
+                                tc_home        = inputs%site_info%tc_home, &
                                 c4             = .false., &
                                 method_optci   = "prentice14", &
                                 method_jmaxlim = "wang17" &
@@ -270,7 +271,7 @@ subroutine gs_leuning( rad_top, rad_net, tl, ea, lai, &
 
     real,    intent(in)    :: rad_top ! PAR dn on top of the canopy, w/m2
     real,    intent(in)    :: rad_net ! PAR net on top of the canopy, w/m2
-    real,    intent(in)    :: tl   ! leaf temperature, degK
+    real,    intent(in)    :: tl   ! leaf temperature, K
     real,    intent(in)    :: ea   ! specific humidity in the canopy air, kg/kg
     real,    intent(in)    :: lai  ! leaf area index
     !real,    intent(in)    :: leaf_age ! age of leaf since budburst (deciduos), days
@@ -519,14 +520,14 @@ function qscomp(T, p) result(qsat)
     !--------Output
     real :: qsat ! Output type: saturated specific humidity, kg/kg
     !--------Inputs
-    real :: T    ! temperature, degK
+    real :: T    ! temperature, K
     real :: p    ! pressure, Pa
     !--------local var
     real :: myesat ! sat. water vapor pressure
-    real :: Temp ! degC
+    real :: Temp ! deg C
 
     ! calculate saturated specific humidity
-    Temp = T - 273.16 ! degC
+    Temp = T - 273.16 ! deg C
     myesat=MIN(610.78*exp(17.27*Temp/(Temp+237.3)), p) ! Pa
     qsat = 0.622*myesat /(p - 0.378*myesat )
 

src/gpp_pmodel.mod.f90

@@ -44,7 +44,7 @@ module md_gpp_pmodel
 
 contains
 
-  subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd)
+  subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd, tc_home)
     !//////////////////////////////////////////////////////////////////
     ! Wrapper function to call to P-model. 
     ! Calculates meteorological conditions with memory based on daily
@@ -63,7 +63,7 @@ subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd)
     type(gridtype)      :: grid
     logical, intent(in) :: init                              ! is true on the very first simulation day (first subroutine call of each gridcell)
     logical, intent(in) :: in_ppfd                           ! whether to use PPFD from forcing or from SPLASH output
-
+    real, intent(in)    :: tc_home                           ! long-term mean max temp of the warmest month (deg C)
     ! local variables
     type(outtype_pmodel) :: out_pmodel              ! list of P-model output variables
     type(climate_type)   :: climate_acclimation     ! list of climate variables to which P-model calculates acclimated traits
@@ -151,6 +151,7 @@ subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd)
                               tc             = temp_memory, &
                               vpd            = vpd_memory, &
                               patm           = patm_memory, &
+                              tc_home        = tc_home, &
                               c4             = params_pft_plant(pft)%c4, &
                               method_optci   = "prentice14", &
                               method_jmaxlim = "wang17" &
@@ -209,8 +210,17 @@ subroutine gpp( tile, tile_fluxes, co2, climate, grid, init, in_ppfd)
       tile_fluxes(lu)%plant(pft)%iwue    = out_pmodel%iwue
 
       ! quantities with instantaneous temperature response
-      tile_fluxes(lu)%plant(pft)%vcmax = calc_ftemp_inst_vcmax( climate%dtemp, climate%dtemp, tcref = 25.0 ) * out_pmodel%vcmax25
-      tile_fluxes(lu)%plant(pft)%jmax  = calc_ftemp_inst_jmax(  climate%dtemp, climate%dtemp, tcref = 25.0 ) * out_pmodel%jmax25
+      tile_fluxes(lu)%plant(pft)%vcmax = calc_ftemp_inst_vcmax( &
+          tc_leaf   = climate%dtemp, &
+          tc_growth = temp_memory, &
+          ! no tc_home needed for calc_ftemp_inst_vcmax
+          tc_ref    = 25.0 ) * out_pmodel%vcmax25
+
+      tile_fluxes(lu)%plant(pft)%jmax = calc_ftemp_inst_jmax( &
+          tc_leaf   = climate%dtemp, &
+          tc_growth = temp_memory, &
+          tc_home   = tc_home, &
+          tc_ref    = 25.0 ) * out_pmodel%jmax25
 
       !----------------------------------------------------------------
       ! Stomatal conductance

src/interface_biosphere_pmodel.mod.f90

@@ -34,6 +34,7 @@ module md_interface_pmodel
     real                                    :: pco2
     type(gridtype)                          :: grid
     real                                    :: whc_prescr
+    real                                    :: tc_home
     type(climate_type), dimension(ndayyear) :: climate
     type(vegcover_type), dimension(ndayyear):: vegcover
     type(outtype_steering)                  :: steering_state

src/interface_in_biosphere_biomee.mod.f90

@@ -20,7 +20,7 @@ module md_interface_in_biomee
 
   !===== Number of parameters
   integer, public, parameter :: nvars_params_siml    = 11
-  integer, public, parameter :: nvars_site_info      = 3
+  integer, public, parameter :: nvars_site_info      = 4
   integer, public, parameter :: nvars_params_tile    = 19
   integer, public, parameter :: nvars_init_soil      = 4
   integer, public, parameter :: nvars_init_cohorts   = 9
@@ -203,6 +203,7 @@ module md_interface_in_biomee
     real :: lon
     real :: lat
     real :: elv ! elevation
+    real :: tc_home
 
   contains
 
@@ -407,9 +408,10 @@ subroutine populate_site_info(self, site_info)
     real(kind=c_double), dimension(nvars_site_info),   intent(in)  :: site_info
 
     ! Site info
-    self%lon = real( site_info(1) )
-    self%lat = real( site_info(2) )
-    self%elv = real( site_info(3) )
+    self%lon     = real( site_info(1) )
+    self%lat     = real( site_info(2) )
+    self%elv     = real( site_info(3) )
+    self%tc_home = real( site_info(4) )
   end subroutine populate_site_info
 
   subroutine populate_spec_data(self, params_species)