diff --git a/data-raw/generate_biomee_driver_data_DBEN.R b/data-raw/generate_biomee_driver_data_DBEN.R new file mode 100644 index 00000000..8121dd39 --- /dev/null +++ b/data-raw/generate_biomee_driver_data_DBEN.R @@ -0,0 +1,258 @@ +#!/usr/bin/env Rscript + +library(tidyverse) +library(rsofun) + +# load script arguments +args <- commandArgs(trailingOnly = TRUE) + +# load data +if (is.na(args[1])){ + load("data-raw/CH-LAE_forcing.rda") +} else { + load(args[1]) +} + +# sitename +sitename <- "CH-Lae" + +# Take only year 2004 to 2014, corresponding to subset of data for site CH-Lae +siteinfo <- tibble( + sitename = "CH-Lae", + lon = 8.365, + lat = 47.47808, + elv = 700, + year_start = 2004, + year_end = 2014, + classid = NA, + c4 = FALSE, + whc = NA, + koeppen_code = NA, + igbp_land_use = "Mixed Forests", + plant_functional_type = "Broadleaf trees") + +siteinfo <- siteinfo %>% + dplyr::mutate(date_start = lubridate::ymd(paste0(year_start, "-01-01"))) %>% + dplyr::mutate(date_end = lubridate::ymd(paste0(year_end, "-12-31"))) + +# load model parameters (valid ones) +params_siml <- tibble( + spinup = TRUE, + spinupyears = 250, + recycle = 1, + firstyeartrend = 2009, + nyeartrend = 1, + outputhourly = TRUE, + outputdaily = TRUE, + do_U_shaped_mortality = TRUE, + update_annualLAImax = TRUE, + do_closedN_run = TRUE, + method_photosynth = "gs_leuning", + method_mortality = "dbh" +) + +params_siml_pmodel <- params_siml +params_siml_pmodel$method_photosynth <- "pmodel" + +params_tile <- tibble( + soiltype = 3, + FLDCAP = 0.4, + WILTPT = 0.05, + K1 = 2.0, + K2 = 0.05, + K_nitrogen = 8.0, + MLmixRatio = 0.8, + etaN = 0.025, + LMAmin = 0.02, + fsc_fine = 1.0, + fsc_wood = 0.0, + GR_factor = 0.33, + l_fract = 0.0, + retransN = 0.0, + f_initialBSW = 0.2, + f_N_add = 0.02, + tf_base = 1, + par_mort = 1, + par_mort_under = 1 +) + +params_species <- tibble( + lifeform = rep(1,16), + phenotype = c(0,1,1,rep(1,13)), + pt = rep(0,16), + alpha_FR = rep(1.2,16), + rho_FR = rep(200,16), + root_r = rep(2.9E-4,16), + root_zeta = rep(0.29,16), + Kw_root = rep(3.5e-09,16), + leaf_size = rep(0.04,16), + Vmax = rep(35.0E-6,16), + Vannual = rep(1.2,16), + wet_leaf_dreg = rep(0.3,16), + m_cond = rep(7.0,16), + alpha_phot = rep(0.06,16), + gamma_L = rep(0.02,16), + gamma_LN = rep(70.5 ,16), + gamma_SW = rep(0.08,16), + gamma_FR = rep(12.0,16), + tc_crit = rep(283.16,16), + tc_crit_on = rep(280.16,16), + gdd_crit = rep(280.0,16), + seedlingsize = rep(0.05,16), + LNbase = rep(0.8E-3,16), + lAImax = rep(3.5,16), + Nfixrate0 = rep(0,16), + NfixCost0 = rep(12,16), + phiCSA = rep(0.25E-4,16), + mortrate_d_c = rep(0.01,16), + mortrate_d_u = rep(0.075,16), + maturalage = rep(5,16), + fNSNmax = rep(5,16), + LMA = c(0.05,0.17,0.11,rep(0.1,13)), + rho_wood = c(590,370,350,rep(300,13)), + alphaBM = rep(5200,16), + thetaBM = c(2.36,2.30,2.54,rep(2.30,13)), + kphio = rep(0.05,16), + phiRL = rep(3.5,16), + LAI_light = rep(3.5,16) +) + +params_soil <- tibble( + type = c("Coarse","Medium","Fine","CM","CF","MF","CMF","Peat","MCM"), + GMD = c(0.7, 0.4, 0.3, 0.1, 0.1, 0.07, 0.007, 0.3, 0.3), + GSD = c(5.0, 5.3, 7.4, 6.1, 6.1, 14.0, 15.0, 7.4, 7.4), + vwc_sat = c(0.380, 0.445, 0.448, 0.412, 0.414, 0.446, 0.424, 0.445, 0.445), + chb = c(3.5,6.4,11.0,4.8,6.3,8.4,6.3,6.4,6.4), + psi_sat_ref = c(-600, -790, -910, -1580, -1680, -1880, -5980, -790, -790), + k_sat_ref = c(130.8, 75.1, 53.2, 12.1, 11.1, 12.7, 1.69, 53.2, 53.2), + alphaSoil = rep(1, 9), + heat_capacity_dry = c(1.2e6, 1.1e6, 1.1e6, 1.1e6, 1.1e6, 1.1e6, 1.1e6, 1.4e6, 1.0) +) + +init_cohort <- tibble( + init_cohort_species = rep(1, 10), + init_cohort_nindivs = rep(0.05,10), + init_cohort_bsw = rep(0.05,10), + init_cohort_bHW = rep(0.0, 10), + init_cohort_nsc = rep(0.05,10) +) + +init_soil <- tibble( #list + init_fast_soil_C = 0.0, + init_slow_soil_C = 0.0, + init_Nmineral = 0.015, + N_input = 0.0008 +) + +#---- gs leuning formatting ----- +forcing <- forcingLAE %>% + dplyr::group_by( + lubridate::month(datehour), + lubridate::day(datehour), + lubridate::hour(datehour)) %>% + summarise_at(vars(1:13), list(~mean(., na.rm = TRUE))) %>% + rename(month=`lubridate::month(datehour)`,day=`lubridate::day(datehour)`) %>% + ungroup() +forcing <- forcing %>% + rename(year=YEAR, doy= DOY, hour=HOUR, par=PAR, ppfd=Swdown, temp=TEMP, temp_soil=SoilT, rh=RH, + prec=RAIN, wind=WIND, patm=PRESSURE, co2=aCO2_AW, swc=SWC) %>% + mutate(snow=NA,vpd=NA, ccov_int=NA,ccov=NA, + date = make_date(year,month,day)) %>% + select(-c(month,day)) %>% + relocate(date,year,doy,hour,temp,temp_soil,prec,snow,vpd,rh,ppfd,par,patm,wind,ccov_int,ccov,co2,swc) +#forcing <- bind_rows(replicate(2, forcing, simplify = FALSE)) + +biomee_gs_leuning_drivers <- tibble( + sitename, + site_info = list(tibble(siteinfo)), + params_siml = list(tibble(params_siml)), + params_tile = list(tibble(params_tile)), + params_species = list(tibble(params_species)), + params_soil = list(tibble(params_soil)), + init_cohort = list(tibble(init_cohort)), + init_soil = list(tibble(init_soil)), + forcing = list(tibble(forcing)) + ) + +save(biomee_gs_leuning_drivers, + file ="data/biomee_gs_leuning_drivers.rda", + compress = "xz") + +#---- p-model formatting ----- +forcing <- forcingLAE %>% + dplyr::group_by( + lubridate::month(datehour), + lubridate::day(datehour)) %>% + summarise_at(vars(1:13), + list(~mean(., na.rm = TRUE))) %>% + rename(month=`lubridate::month(datehour)`,day=`lubridate::day(datehour)`) %>% + ungroup() +forcing <- forcing %>% + rename(year=YEAR, doy= DOY, hour=HOUR, par=PAR, ppfd=Swdown, temp=TEMP, temp_soil=SoilT, rh=RH, + prec=RAIN, wind=WIND, patm=PRESSURE, co2=aCO2_AW, swc=SWC) %>% + mutate(snow=NA,vpd=NA, ccov_int=NA,ccov=NA, + date = make_date(year,month,day)) %>% + select(-c(month,day)) %>% + relocate(date,year,doy,hour,temp,temp_soil,prec,snow,vpd,rh,ppfd,par,patm,wind,ccov_int,ccov,co2,swc) +#forcing <- bind_rows(replicate(2, forcing, simplify = FALSE)) + +biomee_p_model_drivers <- tibble( + sitename, + site_info = list(tibble(siteinfo)), + params_siml = list(tibble(params_siml_pmodel)), + params_tile = list(tibble(params_tile)), + params_species = list(tibble(params_species)), + params_soil = list(tibble(params_soil)), + init_cohort = list(tibble(init_cohort)), + init_soil = list(tibble(init_soil)), + forcing =list(tibble(forcing)) + ) + +save(biomee_p_model_drivers, + file ="data/biomee_p_model_drivers.rda", + compress = "xz") + +# run the model gs-leuning +biomee_gs_leuning_output <- runread_biomee_f( + biomee_gs_leuning_drivers, + makecheck = TRUE, + parallel = FALSE +)$data[[1]]$output_annual_tile[c('year','GPP','plantC')] + +cowplot::plot_grid( + biomee_gs_leuning_output %>% + ggplot() + + geom_line(aes(x = year, y = GPP)) + + theme_classic()+labs(x = "Year", y = "GPP"), + biomee_gs_leuning_output %>% + ggplot() + + geom_line(aes(x = year, y = plantC)) + + theme_classic()+labs(x = "Year", y = "plantC") +) + +save(biomee_gs_leuning_output, + file ="data/biomee_gs_leuning_output.rda", + compress = "xz") + +# run the model p-model +biomee_p_model_output <- runread_biomee_f( + biomee_p_model_drivers, + makecheck = TRUE, + parallel = FALSE +)$data[[1]]$output_annual_tile[c('year','GPP','plantC')] + +cowplot::plot_grid( + biomee_p_model_output %>% + ggplot() + + geom_line(aes(x = year, y = GPP)) + + theme_classic()+labs(x = "Year", y = "GPP"), + biomee_p_model_output %>% + ggplot() + + geom_line(aes(x = year, y = plantC)) + + theme_classic()+labs(x = "Year", y = "plantC") +) + +save(biomee_p_model_output, + file = "data/biomee_p_model_output.rda", + compress = "xz") +