diff --git a/.Rbuildignore b/.Rbuildignore index e6d99f24..2ea4d70d 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -27,3 +27,6 @@ ^.*-logo\..*$ ^codemeta\.json$ ^CRAN-SUBMISSION$ + +# Debugging +^\.vscode$ diff --git a/.github/workflows/test-coverage-debug.yaml b/.github/workflows/test-coverage-debug.yaml new file mode 100644 index 00000000..be54bb65 --- /dev/null +++ b/.github/workflows/test-coverage-debug.yaml @@ -0,0 +1,39 @@ +on: + push: + branches: + - main + - master + - phydro + pull_request: + branches: + - main + - master + - phydro + +name: test-coverage-debug +jobs: + test-coverage-debug: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - uses: actions/checkout@v3 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::devtools + + - name: Install package + run: | + devtools::install() + shell: Rscript {0} + + - name: Run MWE + run: | + source(".vscode/debug.R") + shell: Rscript {0} diff --git a/.vscode/debug.R b/.vscode/debug.R new file mode 100644 index 00000000..3441fab0 --- /dev/null +++ b/.vscode/debug.R @@ -0,0 +1,121 @@ +# NOTE: This R script is called by the debugger (launch.json or through CLI) and +# the function check_NA_output() must be modified as a reproducible example +# showcasing the issue to be debugged + +# ================================================================= ---- +# Further reading: +# Following: https://blog.davisvaughan.com/posts/2022-03-11-using-vs-code-to-debug-r-packages-with-c-code/#setup-function +# which is a detailed description (but only for mac using llvm) +# and: https://www.maths.ed.ac.uk/~swood34/RCdebug/RCdebug.html +# and: https://tdhock.github.io/blog/2019/gdb/ + +# these links: https://github.com/renkun-ken/vscode-rcpp-demo and https://github.com/renkun-ken/vscode-cpp11-demo +# provide an alternative example project (that I have not tested) + +# ================================================================= ---- +# Method 1) Debug this locally in VSCode by launching the debugger through 'launch.json' and set breakpoints +# Method 2) +# To use the exact same configuration as in the GitHub actions use: nektos/act or its implementation in form of +# a VSCode extension 'GitHub Local Actions', making use of Docker containers. +# On the CLI: a) run: ~/bin/act -r -W '.github/workflows/test-coverage-debug.yaml' +# and b) connect: docker exec -it `docker ps -q | head -n1` bash +# and c) debug code in Virtual Machine +# Or (recommended) through VSCode extension: +# a) GitHub Local Actions: make sure to activate the Option reuse=true and run the container +# b) Remote Explorer and Dev Containers: connect VSCode to the running container +# c) debug code in Virtual Machine +# For both approaches step c) requires to install on the virtual machine: +# GDB: sudo apt-get install gdb +# FortLS: pipx install fortls +# and the necessary VSCode extensions incl. path to fortls /root/.local/bin/fortls instead of /home/fabian/.local/bin/fortls). +# ================================================================= ---- + + +# # Below do +# devtools::clean_dll() +# devtools::load_all() +# # to re-build and re-install the R pkg from the working directory. +# # Then modify the code in check_NA_output() to create a reproducibe example showcasing the bug. +# # Launch the debugger with launch.json and set breakpoints. + + +check_NA_output <- function(flag = TRUE){ + if(flag){ + # DOING THIS GIVES NAN: + devtools::clean_dll() |> print() + withr::with_options(list(), + devtools::load_all()) + # likely because it builds rsofun with: + # gfortran -fpic -g -O2 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -ffile-prefix-map=/build/r-base-4sITk6/r-base-4.4.2=. + # -fstack-protector-strong -fstack-clash-protection -fcf-protection -fdebug-prefix-map=/build/r-base-4sITk6/r-base-4.4.2=/usr/src/r-base-4.4.2-1.2404.0 + # -g -O0 -c params_core.mod.f90 -o params_core.mod.o + # NOTE HOW THIS IS *having* -g -O0 on the last line before -c + # see here: https://forum.posit.co/t/getting-rid-of-debug-flags-in-devtools-load-all/186624 + } else { + # DOING THIS GIVES NO NAN: + devtools::clean_dll() |> print() + withr::with_options(list(pkg.build_extra_flags = FALSE), + devtools::load_all()) + # likely because it builds rsofun with: + # gfortran -fpic -g -O2 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -ffile-prefix-map=/build/r-base-4sITk6/r-base-4.4.2=. + # -fstack-protector-strong -fstack-clash-protection -fcf-protection -fdebug-prefix-map=/build/r-base-4sITk6/r-base-4.4.2=/usr/src/r-base-4.4.2-1.2404.0 + # -c params_core.mod.f90 -o params_core.mod.o + # NOTE HOW THIS IS *not having* -g -O0 on the last line before -c + } + + # THIS YIELDS NA UNDER CERTAIN CONDITIONS + params_modl_phydro <- list( + kphio = 0.04998, + kphio_par_a = 0.01, # set to zero to disable temperature-dependence of kphio + kphio_par_b = 1.0, + rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous + tau_acclim = 30.0, + kc_jmax = 0.41, + phydro_K_plant = 5e-17, # TODO: add documentaiton: Phydro: Plant conductivity + phydro_p50_plant = -0.46, # TODO: add documentaiton: Phydro: Plant P50 + phydro_b_plant = 1, # TODO: add documentaiton: Phydro: shape parameter of vulnerability curve + phydro_alpha = 0.08, # TODO: add documentaiton: Phydro: Cost of Jmax + phydro_gamma = 0.065, # TODO: add documentaiton: Phydro: Cost of hydraulics + bsoil = 3, # TODO: add documentaiton: Phydro: parameter converting RZWSC to predawn water potential (depends on rooting system hence PFT specific) + Ssoil = 113 # TODO: add documentaiton: Phydro: parameter converting RZWSC to predawn water potential (depends on rooting system hence PFT specific) + ) + + # read in demo data + df_drivers <- rsofun::p_model_drivers_formatPhydro + + mod4 <- run_pmodel_f_bysite( + sitename = df_drivers$sitename[1], + params_siml = dplyr::mutate(df_drivers$params_siml[[1]], use_phydro = TRUE, use_pml = TRUE, use_gs = TRUE), + site_info = mutate(df_drivers$site_info[[1]], whc = 253), + forcing = df_drivers$forcing[[1]], + forcing_acclim = df_drivers$forcing[[1]], + params_modl = params_modl_phydro, + makecheck = TRUE + ) + # print(df_drivers$params_siml[[1]]) + # print(slice(tibble(mod4), c(1, 70, 1200, 1400, 2000, 2180))) + print(slice(tibble(mod4), c(1, 70, 1200, 1400, 2000, 2180)) |> select(gpp, aet, le, wscal, wcont)) + + # THIS OUTPUTS SOMETIMES: + # gpp aet le wscal wcont + # + # 1 1.37 0.109 269513. NaN NaN + # 2 0.759 1.36 3348090. NaN NaN + # 3 0.917 2.76 6818759 NaN NaN + # 4 0.782 0.326 803802. NaN NaN + # 5 0.363 5.73 13932128 NaN NaN + # 6 1.20 -0.323 -799910 NaN NaN + # BUT OTHER TIMES IT PROVIDES: + # gpp aet le wscal wcont + # + # 1 2.48 0.109 269531. 0.627 159. + # 2 4.29 1.36 3348880. 0.828 210. + # 3 6.38 2.76 6819973 0.964 244. + # 4 1.12 0.326 803872. 0.414 105. + # 5 4.19 5.73 13935005 0.442 112. + # 6 1.46 -0.323 -799917 1 253 +} +# check_NA_output(flag = TRUE) # Returns NaN: check why NaN appear with breakpoints at break gpp_pmodel.mod.f90:297 if count >= 325+1 + # and then see how e (sapflux?) is set to 0, leading to NaN on that day at photosynth_phydro.mod.f90:1544 +check_NA_output(flag = FALSE) # Returns no NaN: check what is computed with breakpoints at break gpp_pmodel.mod.f90:297 if count >= 325+1 + diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 00000000..7d5ebe4c --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,51 @@ +{ + "version": "0.2.0", + "configurations": [ + { + "name": "Debug R-Fortran", + "type": "cppdbg", + "request": "launch", + "program": "/usr/lib/R/bin/exec/R", + "args": [ + "--vanilla", + "-e", + "source('.vscode/debug.R')" + ], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [ + {"name":"R_HOME", "value":"/usr/lib/R"}], + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ], + // "preLaunchTask": "build R pkg (for debugging)" + }, + { + "name": "Docker debug R-Fortran", + "type": "cppdbg", + "request": "launch", + "program": "/opt/R/4.4.3/lib/R/bin/exec/R", + "args": [ + "--vanilla", + "-e", + "source('.vscode/debug.R')" + ], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [ + {"name":"R_HOME", "value":"/opt/R/4.4.3/lib/R"}, + {"name":"LD_LIBRARY_PATH", "value":"/opt/R/4.4.3/lib/R/lib/"}], + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ] + } + ] +} \ No newline at end of file diff --git a/R/calib_sofun.R b/R/calib_sofun.R index ba620806..0da1c567 100644 --- a/R/calib_sofun.R +++ b/R/calib_sofun.R @@ -53,7 +53,7 @@ #' rd_to_vcmax = 0.014, #' tau_acclim = 30, #' kc_jmax = 0.41, -#' whc = 432 +#' gw_calib = 2.0 #' ) #' #' # Define calibration settings diff --git a/R/cost_likelihood_phydro.R b/R/cost_likelihood_phydro.R index d60b2626..10dad334 100644 --- a/R/cost_likelihood_phydro.R +++ b/R/cost_likelihood_phydro.R @@ -71,7 +71,7 @@ #' phydro_alpha = 0.08, #' bsoil = 3, #' Ssoil = 113, -#' whc = 253, +#' gw_calib = 2.0, #' # kphio = 0.09423773, # setup ORG in Stocker et al. 2020 GMD #' # kphio_par_a = 0.0, # set to zero to disable temperature-dependence of kphio #' # kphio_par_b = 1.0, diff --git a/R/cost_likelihood_pmodel.R b/R/cost_likelihood_pmodel.R index 20144f6d..2f07cc7f 100644 --- a/R/cost_likelihood_pmodel.R +++ b/R/cost_likelihood_pmodel.R @@ -71,7 +71,7 @@ #' rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous #' tau_acclim = 30.0, #' kc_jmax = 0.41, -#' whc = 430 +#' gw_calib = 2.0 #' ) #' ) cost_likelihood_pmodel <- function( diff --git a/R/cost_rmse_pmodel.R b/R/cost_rmse_pmodel.R index 93139472..fab7b4c5 100644 --- a/R/cost_rmse_pmodel.R +++ b/R/cost_rmse_pmodel.R @@ -62,7 +62,7 @@ #' rd_to_vcmax = 0.014, # from Atkin et al. 2015 for C3 herbaceous #' tau_acclim = 30.0, #' kc_jmax = 0.41, -#' whc = 240 +#' gw_calib = 2.0 #' ) #' ) diff --git a/R/run_pmodel_f_bysite.R b/R/run_pmodel_f_bysite.R index 448b4b4c..cf23bbca 100644 --- a/R/run_pmodel_f_bysite.R +++ b/R/run_pmodel_f_bysite.R @@ -148,7 +148,7 @@ #' params_siml = p_model_drivers_formatPhydro$params_siml[[1]], #' site_info = p_model_drivers_formatPhydro$site_info[[1]], #' forcing = p_model_drivers_formatPhydro$forcing[[1]], -#' forcing_acclim = p_model_drivers_formatPhydro$forcing_daytime[[1]], +#' forcing_acclim = p_model_drivers_formatPhydro$forcing_daytime[[1]] |> dplyr::mutate(vwind=2.0), # TODO: update p_model_drivers_formatPhydro #' params_modl = params_modl #' ) run_pmodel_f_bysite <- function( # TODO: Above docstring appears duplicated in runread_pmodel_f.R. This redunduncy should be reduced. @@ -421,6 +421,7 @@ run_pmodel_f_bysite <- function( # TODO: Above docstring appears duplicated in r as.numeric(params_modl$rd_to_vcmax), as.numeric(params_modl$tau_acclim), as.numeric(params_modl$kc_jmax), + as.numeric(params_modl$gw_calib), ifelse(params_siml$use_phydro, no = dummy_val, yes = params_modl$phydro_K_plant), @@ -558,6 +559,7 @@ run_pmodel_f_bysite <- function( # TODO: Above docstring appears duplicated in r required_param_names <- list( phydro_model = c( # P-hydro model needs these parameters: 'bsoil', + 'gw_calib', 'kc_jmax', 'kphio', 'kphio_par_a', @@ -573,6 +575,7 @@ required_param_names <- list( ), p_model = c(# P-model needs these parameters: 'beta_unitcostratio', + 'gw_calib', 'kc_jmax', 'kphio', 'kphio_par_a', diff --git a/src/Makevars b/src/Makevars index ce6592f5..d221a0cd 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,71 +1,65 @@ +################################# +# For debugging purpose, the following can be added to your local ~/.R/Makevars file +#PKG_FFLAGS = -frecursive -fbounds-check -fcheck=all -Wall -Wextra -pedantic -g -O0 -fbacktrace -ffpe-trap=invalid,zero,overflow -finit-real=snan -finit-integer=-9999999 -finit-derived +# -frecursive: to avoid "Warning: Array 'out_biosphere' at (1) is larger than limit ..." +# -fbounds-check: https://scicomp.stackexchange.com/questions/36045/ifort-everithing-ok-but-with-gfortran-segmentation-fault +# -Wall -Wextra -Wpedantic: to get more warnings +# -fbacktrace: Specify that, when a runtime error is encountered, the Fortran runtime library should output a backtrace of the error. This option only has influence for compilation of the Fortran main program. +# -fcheck=all: enable "run-time tests", such as, for instance, array bounds checks. +# -ffpe-trap=invalid,zero,overflow: Traps usual FP exceptions (for better debugging) +# The following are useful for debugging uninitialized values. It should not be used in production! +# -finit-real=snan -finit-integer=-9999999: Initialize uninitialized values to NaN for reals, and -9999999 for integers. +# -finit-derived: applies the settings above to fields in derive types too + +# Link-time optimization, see https://cran.r-project.org/doc/manuals/R-admin.html#LTO-with-GCC-1 +#LTO_OPT = -flto +#LTO_FC_OPT = -flto +#IN ADDITION,to adding these flags to ~/.R/Makevars, run 'R CMD INSTALL --use-LTO ' +################################## + # C objects C_OBJS = wrappersc.o -# Fortran objects: refer to file names , order reflects dependency structure -FT_OBJS = params_core.mod.o sofunutils.mod.o grid_siterun.mod.o params_siml_pmodel.mod.o params_siml_biomee.mod.o forcing_siterun_pmodel.mod.o forcing_siterun_biomee.mod.o params_soil_biomee.mod.o \ - interface_biosphere_pmodel.mod.o interface_biosphere_biomee.mod.o tile_pmodel.mod.o plant_pmodel.mod.o soiltemp_sitch.mod.o waterbal_splash.mod.o vegdynamics_pmodel.mod.o \ - photosynth_pmodel.mod.o photosynth_phydro.mod.o \ - gpp_pmodel.mod.o gpp_biomee.mod.o biosphere_pmodel.mod.o biosphere_biomee.mod.o vegetation_biomee.mod.o soil_biomee.mod.o sofun_r.o +# Fortran objects: refer to file names , order reflects dependency structure (Mayeuls solution from main) +FT_OBJS = pmodel.mod.o biomee.mod.o -all: $(SHLIB) clean +all: $(SHLIB) $(SHLIB): $(FT_OBJS) $(C_OBJS) -# Dependency of objects (?) +# Source (object) of Fortran modules # : -sofun_r.o: interface_biosphere_pmodel.mod.o interface_biosphere_biomee.mod.o params_core.mod.o params_core.mod.o params_soil_biomee.mod.o params_siml_pmodel.mod.o params_siml_biomee.mod.o +pmodel.mod.o: interface_biosphere_pmodel.mod.o params_core.mod.o \ + biosphere_pmodel.mod.o params_siml_pmodel.mod.o +biomee.mod.o: interface_biosphere_biomee.mod.o biosphere_biomee.mod.o params_core.mod.o interface_biosphere_pmodel.mod.o: forcing_siterun_pmodel.mod.o params_siml_pmodel.mod.o params_core.mod.o -interface_biosphere_biomee.mod.o: forcing_siterun_biomee.mod.o params_soil_biomee.mod.o params_siml_biomee.mod.o params_core.mod.o +interface_biosphere_biomee.mod.o: forcing_siterun_biomee.mod.o params_siml_biomee.mod.o params_core.mod.o params_soil_biomee.mod.o forcing_siterun_pmodel.mod.o: params_core.mod.o params_siml_pmodel.mod.o grid_siterun.mod.o sofunutils.mod.o forcing_siterun_biomee.mod.o: params_core.mod.o params_siml_biomee.mod.o grid_siterun.mod.o params_soil_biomee.mod.o: params_core.mod.o -tile_pmodel.mod.o: params_core.mod.o interface_biosphere_pmodel.mod.o -waterbal_splash.mod.o: params_core.mod.o tile_pmodel.mod.o plant_pmodel.mod.o sofunutils.mod.o +tile_pmodel.mod.o: params_core.mod.o interface_biosphere_pmodel.mod.o plant_pmodel.mod.o +waterbal_splash.mod.o: params_core.mod.o tile_pmodel.mod.o plant_pmodel.mod.o sofunutils.mod.o gpp_pmodel.mod.o gpp_pmodel.mod.o: params_core.mod.o sofunutils.mod.o plant_pmodel.mod.o tile_pmodel.mod.o interface_biosphere_pmodel.mod.o photosynth_pmodel.mod.o photosynth_phydro.mod.o gpp_biomee.mod.o: datatypes.mod.o soil_biomee.mod.o forcing_siterun_biomee.mod.o photosynth_pmodel.mod.o params_core.mod.o sofunutils.mod.o photosynth_pmodel.mod.o: params_core.mod.o sofunutils.mod.o +photosynth_phydro.mod.o: params_core.mod.o sofunutils.mod.o photosynth_pmodel.mod.o soiltemp_sitch.mod.o: params_core.mod.o sofunutils.mod.o tile_pmodel.mod.o interface_biosphere_pmodel.mod.o plant_pmodel.mod.o: params_core.mod.o sofunutils.mod.o interface_biosphere_pmodel.mod.o vegdynamics_pmodel.mod.o: params_core.mod.o tile_pmodel.mod.o plant_pmodel.mod.o gpp_pmodel.mod.o waterbal_splash.mod.o -biosphere_pmodel.mod.o: params_core.mod.o classdefs.mod.o sofunutils.mod.o plant_pmodel.mod.o waterbal_splash.mod.o gpp_pmodel.mod.o vegdynamics_pmodel.mod.o tile_pmodel.mod.o interface_biosphere_pmodel.mod.o soiltemp_sitch.mod.o vegdynamics_pmodel.mod.o -biosphere_biomee.mod.o: params_core.mod.o interface_biosphere_biomee.mod.o datatypes.mod.o soil_biomee.mod.o vegetation_biomee.mod.o -soil_biomee.mod.o: datatypes.mod.o +biosphere_pmodel.mod.o: params_core.mod.o classdefs.mod.o sofunutils.mod.o plant_pmodel.mod.o waterbal_splash.mod.o \ +gpp_pmodel.mod.o vegdynamics_pmodel.mod.o tile_pmodel.mod.o interface_biosphere_pmodel.mod.o soiltemp_sitch.mod.o vegdynamics_pmodel.mod.o +biosphere_biomee.mod.o: params_core.mod.o interface_biosphere_biomee.mod.o datatypes.mod.o soil_biomee.mod.o vegetation_biomee.mod.o \ +soiltemp_sitch.mod.o sofunutils.mod.o +soil_biomee.mod.o: datatypes.mod.o sofunutils.mod.o vegetation_biomee.mod.o: datatypes.mod.o soil_biomee.mod.o gpp_biomee.mod.o datatypes.mod.o: interface_biosphere_biomee.mod.o params_core.mod.o classdefs.mod.o sofunutils.mod.o: params_core.mod.o -photosynth_phydro.mod.o: sofunutils.mod.o photosynth_pmodel.mod.o - -# Source (object) of Fortran modules -# : -sofun_r_mod.mod: sofun_r.o -md_params_core.mod: params_core.mod.o -md_params_siml_pmodel.mod: params_siml_pmodel.mod.o -md_params_siml_biomee.mod: params_siml_biomee.mod.o -md_forcing_pmodel.mod: forcing_siterun_pmodel.mod.o -md_forcing_biomee.mod: forcing_siterun_biomee.mod.o -md_params_soil_biomee.mod: params_soil_biomee.mod.o -md_interface_pmodel.mod: interface_biosphere_pmodel.mod.o -md_interface_biomee.mod: interface_biosphere_biomee.mod.o -md_grid.mod: grid_siterun.mod.o -md_biosphere_pmodel.mod: biosphere_pmodel.mod.o -md_biosphere_biomee.mod: biosphere_biomee.mod.o -md_classdefs.mod: classdefs.mod.o -md_plant_pmodel.mod: plant_pmodel.mod.o -md_waterbal.mod: waterbal_splash.mod.o -md_sofunutils.mod: sofunutils.mod.o -md_tile_pmodel.mod: tile_pmodel.mod.o -md_gpp_pmodel.mod: gpp_pmodel.mod.o -md_gpp_biomee.mod: gpp_biomee.mod.o -md_photosynth.mod: photosynth_pmodel.mod.o -md_photosynth_phydro.mod: photosynth_phydro.mod.o -md_soiltemp.mod: soiltemp_sitch.mod.o -md_vegdynamics_pmodel.mod: vegdynamics_pmodel.mod.o -datatypes.mod: datatypes.o -md_soil_biomee.mod: soil_biomee.o -md_vegetation_biomee.mod: vegetation_biomee.o +params_siml_biomee.mod.o: params_core.mod.o +params_siml_pmodel.mod.o: params_core.mod.o +classdefs.mod.o: params_core.mod.o # Dependency of the C wrapper -wrappersc.o: sofun_r_mod.mod +wrappersc.o: pmodel.mod.o biomee.mod.o clean: - @rm -rf *.mod *.o + @rm -rf *.o *.mod \ No newline at end of file diff --git a/src/gpp_pmodel.mod.f90 b/src/gpp_pmodel.mod.f90 index 7a9c81ee..052eb959 100644 --- a/src/gpp_pmodel.mod.f90 +++ b/src/gpp_pmodel.mod.f90 @@ -17,7 +17,7 @@ module md_gpp_pmodel implicit none private - public params_pft_gpp, gpp, getpar_modl_gpp + public params_pft_gpp, gpp, getpar_modl_gpp, paramstype_gpp, params_gpp !----------------------------------------------------------------------- ! Uncertain (unknown) parameters. Runtime read-in @@ -30,6 +30,7 @@ module md_gpp_pmodel real :: tau_acclim_tempstress real :: par_shape_tempstress real :: kc_jmax + real :: gw_calib end type paramstype_gpp ! PFT-DEPENDENT PARAMETERS @@ -103,7 +104,7 @@ subroutine gpp( tile, tile_fluxes, co2, climate, climate_acclimation, grid, init ! xxx test real :: a_c, a_j, a_returned, fact_jmaxlim - integer, save :: count + integer, save :: count ! TODO: THIS APPEARS NOT TO BE NEEDED. (But since it corresponds to doy, it is useful for conditional breakpoints to debug the NaN issue.) !---------------------------------------------------------------- ! Convert daily mean environmental conditions to conditions to @@ -674,12 +675,15 @@ subroutine getpar_modl_gpp() ! Jmax cost coefficient, c* in Stocker et al., 2020 GMD (Eq 15) and Wang et al., 2017 params_gpp%kc_jmax = myinterface%params_calib%kc_jmax ! 0.41 - + ! Acclimation time scale for photosynthesis (d), multiple lines of evidence suggest about monthly is alright params_gpp%tau_acclim = myinterface%params_calib%tau_acclim ! 30.0 ! Re-interpreted soil moisture stress parameter, previously thetastar = 0.6 params_gpp%soilm_thetastar = myinterface%params_calib%soilm_thetastar + + ! Scaling factor from leaf-level to canopy-level conductance + params_gpp%gw_calib = myinterface%params_calib%gw_calib ! quantum yield efficiency at optimal temperature, phi_0 (Stocker et al., 2020 GMD Eq. 10) params_pft_gpp(:)%kphio = myinterface%params_calib%kphio diff --git a/src/interface_biosphere_pmodel.mod.f90 b/src/interface_biosphere_pmodel.mod.f90 index 1a5de6bb..869391c8 100644 --- a/src/interface_biosphere_pmodel.mod.f90 +++ b/src/interface_biosphere_pmodel.mod.f90 @@ -25,6 +25,7 @@ module md_interface_pmodel real :: rd_to_vcmax real :: tau_acclim real :: kc_jmax + real :: gw_calib real :: phydro_K_plant real :: phydro_p50_plant real :: phydro_b_plant diff --git a/src/sofun_r.f90 b/src/sofun_r.f90 index 9a4bfda3..33f237b6 100644 --- a/src/sofun_r.f90 +++ b/src/sofun_r.f90 @@ -84,7 +84,7 @@ subroutine pmodel_f( & real(kind=c_double), intent(in) :: canopy_height real(kind=c_double), intent(in) :: reference_height integer(kind=c_int), intent(in) :: nt ! number of time steps - real(kind=c_double), dimension(15), intent(in) :: par ! free (calibratable) model parameters + real(kind=c_double), dimension(16), intent(in) :: par ! free (calibratable) model parameters real(kind=c_double), dimension(nt,13), intent(in) :: forcing ! array containing all temporally varying forcing data for instantaneous model (rows: time steps; columns: 1=air temperature, 2=rainfall, 3=vpd, 4=ppfd, 5=net radiation, 6=sunshine fraction, 7=snowfall, 8=co2, 9=fapar, 10=patm, 11=tmin, 12=tmax) real(kind=c_double), dimension(nt,12), intent(in) :: forcing_acclim ! array containing all temporally varying forcing data for acclimating model (rows: time steps; columns: 1=air temperature, 2=rainfall, 3=vpd, 4=ppfd, 5=net radiation, 6=sunshine fraction, 7=snowfall, 8=co2, 9=fapar, 10=patm, 11=tmin, 12=tmax) real(kind=c_double), dimension(nt,23), intent(out) :: output @@ -166,13 +166,14 @@ subroutine pmodel_f( & myinterface%params_calib%rd_to_vcmax = real(par(6)) myinterface%params_calib%tau_acclim = real(par(7)) myinterface%params_calib%kc_jmax = real(par(8)) - myinterface%params_calib%phydro_K_plant = real(par(9)) - myinterface%params_calib%phydro_p50_plant = real(par(10)) - myinterface%params_calib%phydro_b_plant = real(par(11)) - myinterface%params_calib%phydro_alpha = real(par(12)) - myinterface%params_calib%phydro_gamma = real(par(13)) - myinterface%params_calib%bsoil = real(par(14)) - myinterface%params_calib%Ssoil = real(par(15)) + myinterface%params_calib%gw_calib = real(par(9)) + myinterface%params_calib%phydro_K_plant = real(par(10)) + myinterface%params_calib%phydro_p50_plant = real(par(11)) + myinterface%params_calib%phydro_b_plant = real(par(12)) + myinterface%params_calib%phydro_alpha = real(par(13)) + myinterface%params_calib%phydro_gamma = real(par(14)) + myinterface%params_calib%bsoil = real(par(15)) + myinterface%params_calib%Ssoil = real(par(16)) !---------------------------------------------------------------- ! GET VEGETATION COVER (fractional projective cover by PFT) diff --git a/src/tile_pmodel.mod.f90 b/src/tile_pmodel.mod.f90 index 093974a8..091506bd 100644 --- a/src/tile_pmodel.mod.f90 +++ b/src/tile_pmodel.mod.f90 @@ -230,7 +230,7 @@ subroutine initglobal_soil( soil ) ! argument type(soil_type), intent(inout) :: soil - call initglobal_soil_phy( soil%phy ) + call initglobal_soil_phy( soil%phy ) ! TODO: doesn't this to be updated daily, since it represents physical water status of soil... call initglobal_soil_params( soil%params ) soil%phy%wscal = soil%phy%wcont / soil%params%whc @@ -443,15 +443,33 @@ subroutine initdaily_tile_fluxes( tile_fluxes ) tile_fluxes(:)%canopy%drn = 0.0 tile_fluxes(:)%canopy%drnn = 0.0 tile_fluxes(:)%canopy%rnl = 0.0 - tile_fluxes(:)%canopy%dcn = 0.0 + tile_fluxes(:)%canopy%dcn = 0.0 + tile_fluxes(:)%canopy%deet = 0.0 + tile_fluxes(:)%canopy%dpet = 0.0 + tile_fluxes(:)%canopy%dpet_e = 0.0 + ! TODO: missing here, i.e. still nan initialized: deet, dpet, dpet_e tile_fluxes(:)%canopy%daet = 0.0 tile_fluxes(:)%canopy%daet_e = 0.0 tile_fluxes(:)%canopy%daet_soil = 0.0 tile_fluxes(:)%canopy%daet_e_soil = 0.0 tile_fluxes(:)%canopy%daet_canop = 0.0 tile_fluxes(:)%canopy%daet_e_canop = 0.0 + ! TODO: missing here, i.e. still nan initialized: cpa, dtransp, dgs, dgc, + tile_fluxes(:)%canopy%cpa = 0.0 + tile_fluxes(:)%canopy%dtransp = 0.0 + tile_fluxes(:)%canopy%dgs = 0.0 + tile_fluxes(:)%canopy%dgc = 0.0 tile_fluxes(:)%canopy%dgpp = 0.0 tile_fluxes(:)%canopy%drd = 0.0 + ! TODO: missing here, i.e. still nan initialized: assim,vcmax25,jmax25,vcmax,jmax,gs_accl,chi,iwue, + tile_fluxes(:)%canopy%assim = 0.0 + tile_fluxes(:)%canopy%vcmax25 = 0.0 + tile_fluxes(:)%canopy%jmax25 = 0.0 + tile_fluxes(:)%canopy%vcmax = 0.0 + tile_fluxes(:)%canopy%jmax = 0.0 + tile_fluxes(:)%canopy%gs_accl = 0.0 + tile_fluxes(:)%canopy%chi = 0.0 + tile_fluxes(:)%canopy%iwue = 0.0 tile_fluxes(:)%canopy%ppfd_splash = 0.0 tile_fluxes(:)%canopy%dra = 0.0 ! tile_fluxes(:)%canopy%nu = 0.0 @@ -460,6 +478,15 @@ subroutine initdaily_tile_fluxes( tile_fluxes ) do pft = 1,npft tile_fluxes(:)%plant(npft)%dgpp = 0.0 tile_fluxes(:)%plant(npft)%drd = 0.0 + ! TODO: missing here, i.e. still nan initialized: assim,vcmax25,jmax25,vcmax,jmax,gs_accl,chi,iwue, + tile_fluxes(:)%plant(npft)%assim = 0.0 + tile_fluxes(:)%plant(npft)%vcmax25 = 0.0 + tile_fluxes(:)%plant(npft)%jmax25 = 0.0 + tile_fluxes(:)%plant(npft)%vcmax = 0.0 + tile_fluxes(:)%plant(npft)%jmax = 0.0 + tile_fluxes(:)%plant(npft)%gs_accl = 0.0 + tile_fluxes(:)%plant(npft)%chi = 0.0 + tile_fluxes(:)%plant(npft)%iwue = 0.0 tile_fluxes(:)%plant(npft)%dtransp = 0.0 tile_fluxes(:)%plant(npft)%dlatenth = 0.0 tile_fluxes(:)%plant(npft)%dpsi = 0.0 @@ -542,6 +569,7 @@ subroutine diag_daily( tile, tile_fluxes ) ! Sum over PFTs to get canopy-level quantities !---------------------------------------------------------------- do lu=1,nlu + ! TODO: missing here, DEBUG FABIAN: there are many more members in tile_fluxes(lu)%canopy than what is aggregated below. Is this right? tile_fluxes(lu)%canopy%dgpp = sum(tile_fluxes(lu)%plant(:)%dgpp * tile(lu)%plant(:)%fpc_grid) tile_fluxes(lu)%canopy%dtransp = sum(tile_fluxes(lu)%plant(:)%dtransp * tile(lu)%plant(:)%fpc_grid) tile_fluxes(lu)%canopy%drd = sum(tile_fluxes(lu)%plant(:)%drd * tile(lu)%plant(:)%fpc_grid) diff --git a/src/waterbal_splash.mod.f90 b/src/waterbal_splash.mod.f90 index 782a2946..35f3c0f0 100644 --- a/src/waterbal_splash.mod.f90 +++ b/src/waterbal_splash.mod.f90 @@ -11,6 +11,7 @@ module md_waterbal use md_grid, only: gridtype use md_interface_pmodel, only: myinterface use md_sofunutils, only: radians, dgsin, dgcos, degrees + use md_gpp_pmodel, only: params_gpp implicit none @@ -87,7 +88,7 @@ subroutine waterbal( tile, tile_fluxes, grid, climate, fapar, using_phydro, usin !--------------------------------------------------------- ! Canopy transpiration and soil evaporation !--------------------------------------------------------- - call calc_et( tile_fluxes(lu), grid, climate, sw, fapar(lu), using_phydro, using_gs, using_pml ) + call calc_et( tile_fluxes(lu), grid, climate, sw, fapar(lu), using_phydro, using_gs, using_pml, params_gpp%gw_calib ) !--------------------------------------------------------- ! Update soil moisture and snow pack @@ -296,7 +297,7 @@ subroutine solar( tile_fluxes, grid, climate, doy, in_netrad ) end subroutine solar - subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_gs, using_pml ) + subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_gs, using_pml, gw_calib ) !///////////////////////////////////////////////////////////////////////// ! !------------------------------------------------------------------------- @@ -308,7 +309,8 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g type(gridtype), intent(in) :: grid type(climate_type), intent(in) :: climate real, intent(in) :: sw ! evaporative supply rate, mm/hr - real, intent(in) :: fapar + real, intent(in) :: fapar + real, intent(in) :: gw_calib ! scaling factor from leaf-level to canopy-level conductance logical, intent(in) :: using_phydro logical, intent(in) :: using_gs ! Should Pmodel/Phydro gs be used in ET calc? (otherwise, PT formulation will be used) logical, intent(in) :: using_pml ! If using Pmodel/Phydro gs, should ET be calculated using PM equation (otherwise, diffusion equation will be used) @@ -481,12 +483,18 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g ! print*,'in waterbal: gs_accl ', tile_fluxes%canopy%gs_accl ! introduced scaling with LAI lai_fapar = -1.0 / k_beer * log(1.0 - fapar) - gw = max(2.0 * lai_fapar * tile_fluxes%canopy%gs_accl * 1.6 * kR * (climate%dtemp + kTkelvin), eps) + gw = max(gw_calib * lai_fapar * tile_fluxes%canopy%gs_accl * 1.6 * kR * (climate%dtemp + kTkelvin), eps) ! latent energy flux from canopy (W m-2) ! See also calc_transpiration_pm() in photosynth_phydro.mod.f90 tile_fluxes%canopy%daet_e_canop = (epsilon * fapar * tile_fluxes%canopy%drn + (rho_water * cp / gamma) & * ga * climate%dvpd) / (epsilon + 1.0 + ga / gw) + + + ! canopy conductance assuming gw = infinite + tile_fluxes%canopy%dpet_e =(epsilon * fapar * tile_fluxes%canopy%drn + (rho_water * cp / gamma) & + * ga * climate%dvpd) / (epsilon + 1.0) + tile_fluxes%canopy%dpet = dpet_soil + tile_fluxes%canopy%dpet_e * energy_to_mm ! print*,'-----------------------' ! print*,'canopy_height ', myinterface%canopy_height diff --git a/tests/testthat/test-calibration-pmodel.R b/tests/testthat/test-calibration-pmodel.R index 6e276e7b..295744f8 100644 --- a/tests/testthat/test-calibration-pmodel.R +++ b/tests/testthat/test-calibration-pmodel.R @@ -15,7 +15,8 @@ test_that("test GPP calibration routine p-model (BT, likelihood maximization)", beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) settings <- list( @@ -86,7 +87,7 @@ test_that("test GPP calibration routine p-model (GenSA, rmse, all params)", { settings = settings, optim_out = FALSE, # extra arguments for the cost function - par_fixed = list(), + par_fixed = list(gw_calib = 2.0), targets = 'gpp' ) @@ -143,7 +144,7 @@ test_that("test Vcmax25 calibration routine p-model (BT, likelihood, all params) settings = settings, optim_out = FALSE, # arguments for cost function - par_fixed = list(), + par_fixed = list(gw_calib = 2.0), targets = 'vcmax25' ) # plot(pars$mod) @@ -179,7 +180,8 @@ test_that("test Vcmax25 calibration routine p-model (GenSA, rmse)", { beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous # tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) settings <- list( @@ -240,7 +242,8 @@ test_that("test joint calibration routine p-model (BT, likelihood maximization)" beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) settings <- list( diff --git a/tests/testthat/test-model-runs.R b/tests/testthat/test-model-runs.R index 561953b7..95c435ed 100644 --- a/tests/testthat/test-model-runs.R +++ b/tests/testthat/test-model-runs.R @@ -11,7 +11,7 @@ test_that("run_pmodel_f_bysite()", { rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41, - + gw_calib = 2.0, soilm_thetastar = 0.6 * 240, # to recover old setup with soil moisture stress beta_unitcostratio = 146.0 ) @@ -22,7 +22,7 @@ test_that("run_pmodel_f_bysite()", { rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41, - + gw_calib = 2.0, phydro_K_plant = 5e-17, # TODO: add documentaiton: Phydro: Plant conductivity phydro_p50_plant = -0.46, # TODO: add documentaiton: Phydro: Plant P50 phydro_b_plant = 1, # TODO: add documentaiton: Phydro: shape parameter of vulnerability curve @@ -203,7 +203,7 @@ test_that("run_pmodel_f_bysite()", { gpp = c(1.65618813037872, 6.02679443359375, 6.72385692596436, 1.84405922889709, 9.40890026092529, 0.896598398685455), aet = c(0.0496958047151566, 0.866879522800446, 1.87889838218689, 0.156055122613907, 4.60135173797607, -0.136744096875191), le = c(122955.984375, 2140935.5, 4637295, 384523.84375, 11185824, -339057.375), - pet = c(0.103817254304886, 1.45484209060669, 3.08336925506592, 0.335356384515762, 6.51526165008545, -0.377504140138626), + pet = c(0.108938276767731, 1.35599589347839, 2.76327228546143, 0.326252281665802, 5.732346534729, -0.322612106800079), vcmax = c(9.54577535594581e-06, 1.18804200610612e-05, 1.91590133908903e-05, 1.37124443426728e-05, 5.95575438637752e-05, 5.93948425375856e-06), jmax = c(3.08439557556994e-05, 3.66509229934309e-05, 5.80160958634224e-05, 3.691642996273e-05, 0.000109297579911072, 2.01222374016652e-05), vcmax25 = c(3.51696653524414e-05, 3.77493124688044e-05, 5.76805359742139e-05, 3.59945297532249e-05, 5.17318439960945e-05, 2.60039796557976e-05), @@ -231,7 +231,7 @@ test_that("run_pmodel_f_bysite()", { gpp = c(2.47916746139526, 4.29396057128906, 6.37855243682861, 1.12339150905609, 4.18545007705688, 1.46236681938171), aet = c(0.108937680721283, 1.35598421096802, 2.76325654983521, 0.32624351978302, 5.73224258422852, -0.3226118683815), le = c(269530.59375, 3348879.5, 6819973, 803872.375, 13935005, -799917), - pet = c(0.103817254304886, 1.45484209060669, 3.08336925506592, 0.335356384515762, 6.51526165008545, -0.377504140138626), + pet = c(0.108938276767731, 1.35599589347839, 2.76327228546143, 0.326252281665802, 5.732346534729, -0.322612106800079), vcmax = c(6.111431048339e-06, 7.75509124650853e-06, 1.19958594950731e-05, 3.18529146170476e-06, 4.026747046737e-05, 4.28107478001039e-06), jmax = c(2.28475473704748e-05, 2.86447157122893e-05, 4.33829118264839e-05, 5.85015914111864e-06, 7.40396353648975e-05, 1.724714456941e-05), vcmax25 = c(1.77550300577423e-05, 2.02352257474558e-05, 2.96198322757846e-05, 7.76460092311027e-06, 3.53926807292737e-05, 1.51183612615569e-05), @@ -250,7 +250,6 @@ test_that("run_pmodel_f_bysite()", { le_soil = c(98347.3515625, 1303413.875, 2931274, 273886.84375, 5190479, -290766.28125), dpsi = c(0.338661968708038, 1.0493232011795, 1.22259771823883, 0.456368923187256, 1.90744316577911, 0.0780341103672981), psi_leaf = c(-0.624942302703857, -1.11408090591431, -1.22973167896271, -1.60701656341553, -2.70482659339905, -0.0780341103672981)) - expect_equal(dplyr::slice(tibble(mod1), c(1, 70, 1200, 1400, 2000, 2180)), ref1, tolerance = 1e-6) expect_equal(dplyr::slice(tibble(mod2), c(1, 70, 1200, 1400, 2000, 2180)), ref2, tolerance = 1e-6) expect_equal(dplyr::slice(tibble(mod3), c(1, 70, 1200, 1400, 2000, 2180)), ref3, tolerance = 1e-6) @@ -274,7 +273,8 @@ test_that("runread_pmodel_f()", { beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) # read in demo data @@ -318,7 +318,8 @@ test_that("p-model run check Vcmax25", { beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) # read in demo data @@ -335,7 +336,6 @@ test_that("p-model run check Vcmax25", { mutate(x, canopy_height = 5, reference_height = 10))) - for (it in seq_along((df_drivers$forcing))){df_drivers$forcing[[it]]$vwind = 2.0} # TODO: add vwind, TODO: this would need to be changed in rsofun example data # run the SOFUN Fortran P-model mod <- rsofun::run_pmodel_f_bysite( @@ -345,7 +345,7 @@ test_that("p-model run check Vcmax25", { forcing = df_drivers$forcing[[1]], forcing_acclim = df_drivers$forcing[[1]], params_modl = params_modl, - makecheck = FALSE + makecheck = TRUE ) # test if the returned values @@ -388,7 +388,8 @@ test_that("phydro-model run check LE and AET", { beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) # read in demo data diff --git a/tests/testthat/test-quantitative-validation.R b/tests/testthat/test-quantitative-validation.R index 7fbcb672..8533429c 100644 --- a/tests/testthat/test-quantitative-validation.R +++ b/tests/testthat/test-quantitative-validation.R @@ -17,7 +17,8 @@ test_that("p-model quantitative check versus observations (FR-Pue)", { beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) df_drivers <- rsofun::p_model_drivers_formatPhydro # TODO: NOT YET UPDATED FOR PHYDRO (still add default phydro_* parameters) # |> diff --git a/vignettes/files/sensitivity_analysis.Rmd__morrisOut.rds b/vignettes/files/sensitivity_analysis.Rmd__morrisOut.rds index 795922a1..c37d9de4 100644 Binary files a/vignettes/files/sensitivity_analysis.Rmd__morrisOut.rds and b/vignettes/files/sensitivity_analysis.Rmd__morrisOut.rds differ diff --git a/vignettes/files/sensitivity_analysis.Rmd__par_calib.rds b/vignettes/files/sensitivity_analysis.Rmd__par_calib.rds index fd273a51..94392972 100644 Binary files a/vignettes/files/sensitivity_analysis.Rmd__par_calib.rds and b/vignettes/files/sensitivity_analysis.Rmd__par_calib.rds differ diff --git a/vignettes/pmodel_use.Rmd b/vignettes/pmodel_use.Rmd index d9d933d7..cca397bf 100644 --- a/vignettes/pmodel_use.Rmd +++ b/vignettes/pmodel_use.Rmd @@ -184,7 +184,8 @@ params_modl <- list( beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) drivers_tmp <- p_model_drivers |> @@ -271,7 +272,8 @@ pars <- calib_sofun( beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) ) ``` @@ -288,7 +290,8 @@ params_modl <- list( beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) # root zone total water holding capacity (whc) is calibrated, but specified as diff --git a/vignettes/pmodel_use_newdata.Rmd b/vignettes/pmodel_use_newdata.Rmd index 512d9714..c6228d5d 100644 --- a/vignettes/pmodel_use_newdata.Rmd +++ b/vignettes/pmodel_use_newdata.Rmd @@ -268,7 +268,8 @@ params_modl <- list( beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, - kc_jmax = 0.41 + kc_jmax = 0.41, + gw_calib = 2.0 ) # run the model for these parameters @@ -622,7 +623,8 @@ params_modl <- list( phydro_b_plant = 1, phydro_alpha = 0.08, bsoil = 3, - Ssoil = 113 + Ssoil = 113, + gw_calib = 2.0 ) # run the model for these parameters diff --git a/vignettes/sensitivity_analysis.Rmd b/vignettes/sensitivity_analysis.Rmd index e9121757..d3e7b202 100644 --- a/vignettes/sensitivity_analysis.Rmd +++ b/vignettes/sensitivity_analysis.Rmd @@ -62,7 +62,8 @@ ll_pmodel <- function( as.list(par_v), # must be a named list obs = rsofun::p_model_validation, # example data from package drivers = rsofun::p_model_drivers_formatPhydro, #TODO rsofun::p_model_drivers is NOT YET UPDATED FOR PHYDRO (a newformat, b add phydro_ parameters) - targets = "gpp" + targets = "gpp", + par_fixed = list() ) } @@ -78,6 +79,7 @@ ll_pmodel( rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous tau_acclim = 30.0, kc_jmax = 0.41, + gw_calib = 2.0, err_gpp = 0.9 # value from previous simulations )) ``` @@ -97,6 +99,7 @@ par_cal_best <- c( rd_to_vcmax = 0.014, tau_acclim = 30.0, kc_jmax = 0.41, + gw_calib = 2.0, err_gpp = 1 ) @@ -110,6 +113,7 @@ par_cal_min <- c( rd_to_vcmax = 0.01, tau_acclim = 7.0, kc_jmax = 0.2, + gw_calib = 1.5, err_gpp = 0.01 ) @@ -123,6 +127,7 @@ par_cal_max <- c( rd_to_vcmax = 0.1, tau_acclim = 60.0, kc_jmax = 0.8, + gw_calib = 2.5, err_gpp = 4 ) ``` @@ -277,7 +282,8 @@ par_calib <- calib_sofun( beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, tau_acclim = 30.0, - kc_jmax = 0.41), + kc_jmax = 0.41, + gw_calib = 2.0), targets = "gpp" ) @@ -316,7 +322,8 @@ par_calib <- calib_sofun( kphio_par_b = 20, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, - tau_acclim = 30.0), + tau_acclim = 30.0, + gw_calib = 2.0), targets = "gpp" ) ``` @@ -433,7 +440,8 @@ run_pmodel <- function(par){ kphio_par_b = 20, beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, - tau_acclim = 30.0) + tau_acclim = 30.0, + gw_calib = 2.0) ) return(out) }