From 5de81fb6d7ed31ef130a433a40e7df64bae69b2f Mon Sep 17 00:00:00 2001 From: Beni Stocker Date: Thu, 15 Aug 2024 11:21:23 +0200 Subject: [PATCH] tested: identical results for gs-coupled ET with diffusion (not PML). --- src/forcing_siterun_pmodel.mod.f90 | 32 ++++++++-- src/gpp_pmodel.mod.f90 | 87 +++++++++++++-------------- src/plant_pmodel.mod.f90 | 2 +- src/tile_pmodel.mod.f90 | 3 - src/waterbal_splash.mod.f90 | 22 +++---- vignettes/pmodel_use_newdata.Rmd | 96 +++++++++++++++--------------- 6 files changed, 130 insertions(+), 112 deletions(-) diff --git a/src/forcing_siterun_pmodel.mod.f90 b/src/forcing_siterun_pmodel.mod.f90 index 8d637378..b41516e8 100644 --- a/src/forcing_siterun_pmodel.mod.f90 +++ b/src/forcing_siterun_pmodel.mod.f90 @@ -185,23 +185,47 @@ function get_fpc_grid( params_siml ) result( fpc_grid_field ) ! get binary information of PFT presence from simulation parameters fpc_grid_field(:) = 0.0 - ! Code below must follow the same structure as in 'plant_pmodel.mod.f90' + ! Code below must follow the same structure as in plant_pmodel.mod.f90, getpar_modl_plant() pft = 0 if ( params_siml%ltre ) then - ! xxx dirty: call all non-grass vegetation types 'TrE', see indeces above + ! Consider all non-grass vegetation types 'TrE', see indeces above pft = pft + 1 + fpc_grid_field(:) = 0.0 + fpc_grid_field(pft) = 1.0 + end if + + if ( params_siml%ltne ) then + ! Consider all non-grass vegetation types 'TNE', see indeces above + pft = pft + 1 + fpc_grid_field(:) = 0.0 + fpc_grid_field(pft) = 1.0 + end if + + if ( params_siml%ltrd ) then + ! Consider all non-grass vegetation types 'TrE', see indeces above + pft = pft + 1 + fpc_grid_field(:) = 0.0 + fpc_grid_field(pft) = 1.0 + end if + + if ( params_siml%ltnd ) then + ! Consider all non-grass vegetation types 'TrE', see indeces above + pft = pft + 1 + fpc_grid_field(:) = 0.0 fpc_grid_field(pft) = 1.0 end if if ( params_siml%lgr3 ) then - ! xxx dirty: call all grass vegetation types 'Gr3' + ! Consider all grass vegetation types 'Gr3' pft = pft + 1 + fpc_grid_field(:) = 0.0 fpc_grid_field(pft) = 1.0 end if if ( params_siml%lgr4 ) then - ! xxx dirty: call all grass vegetation types 'Gr3' + ! Consider all grass vegetation types 'Gr3' pft = pft + 1 + fpc_grid_field(:) = 0.0 fpc_grid_field(pft) = 1.0 end if diff --git a/src/gpp_pmodel.mod.f90 b/src/gpp_pmodel.mod.f90 index 04aa6caa..36b91562 100644 --- a/src/gpp_pmodel.mod.f90 +++ b/src/gpp_pmodel.mod.f90 @@ -212,6 +212,18 @@ subroutine gpp( tile, tile_fluxes, co2, climate, climate_acclimation, grid, init method_optci = "prentice14", & method_jmaxlim = "wang17" & ) + + ! print*,'kphio', kphio_temp + ! print*,'beta', params_gpp%beta + ! print*,'kc_jmax', params_gpp%kc_jmax + ! print*,'ppfd', ppfd_memory + ! print*,'co2', co2_memory + ! print*,'tc', temp_memory + ! print*,'vpd', vpd_memory + ! print*,'patm', patm_memory + ! print*,'c4', params_pft_plant(pft)%c4 + ! print*,'-------------------------------' + else par_cost = par_cost_type(tile(lu)%plant(pft)%phydro_alpha, & tile(lu)%plant(pft)%phydro_gamma) @@ -240,8 +252,6 @@ subroutine gpp( tile, tile_fluxes, co2, climate, climate_acclimation, grid, init par_control = options & ) - ! print *, temp_memory, ppfd_memory*1e6, kphio_temp, vpd_memory - ! print *, out_phydro_acclim%a, out_phydro_acclim%gs, out_phydro_acclim%dpsi end if else @@ -346,17 +356,14 @@ subroutine gpp( tile, tile_fluxes, co2, climate, climate_acclimation, grid, init !---------------------------------------------------------------- ! Stomatal conductance to CO2 !---------------------------------------------------------------- - ! Apply soil moisture stress function to stomatal conductance - tile_fluxes(lu)%plant(pft)%gs_accl = out_pmodel%gs_setpoint * soilmstress - - ! if (.not. use_phydro) then - ! ! Jaideep NOTE: I have applied the soilmstress factor to gs here because it is needed in calculating canopy transpiration - ! tile_fluxes(lu)%plant(pft)%gs_accl = out_pmodel%gs_setpoint * soilmstress - ! else - ! ! Jaideep NOTE: unit of gs_accl here is mol m-2 s-1. - ! ! Jaideep FIXME: It's too complicated to convert it to unit as in pmodel, but should be done at some point - ! tile_fluxes(lu)%plant(pft)%gs_accl = out_phydro_inst%gs - ! end if + if (.not. use_phydro) then + ! Jaideep NOTE: I have applied the soilmstress factor to gs here because it is needed in calculating canopy transpiration + tile_fluxes(lu)%plant(pft)%gs_accl = out_pmodel%gs_setpoint * soilmstress + else + ! Jaideep NOTE: unit of gs_accl here is mol m-2 s-1. + ! Jaideep FIXME: It's too complicated to convert it to unit as in pmodel, but should be done at some point + tile_fluxes(lu)%plant(pft)%gs_accl = out_phydro_inst%gs + end if !---------------------------------------------------------------- ! Water potentials @@ -375,39 +382,27 @@ subroutine gpp( tile, tile_fluxes, co2, climate, climate_acclimation, grid, init ! Density of water, kg/m^3 rho_water = calc_density_h2o( climate%dtemp, climate%dpatm ) - if (use_pml) then - ! We plug Pmodel/Phydro-derived gs into the PM equation to calculate T (note this is uncoupled PM-transpiration) - ! use PFT-specific gs for this calculation: tile_fluxes(lu)%plant(pft)%gs_accl - - ! TODO: Fill this in - ! tile_fluxes(lu)%plant(pft)%dtransp = PM_EQUATION(tile_fluxes(lu)%plant(pft)%gs_accl) - - else - ! We plug Pmodel/Phydro-derived gs into T = 1.6gsD - if (.not. use_phydro) then - ! Using P-model gs - ! Note here that stomatal conductance is already normalized by patm (=gs/patm) so E = 1.6 * (gs/patm) * vpd, which is the same as 1.6 gs (vpd/patm) - ! but it is expressed per unit absorbed light, so multiply by PPFD*fapar - ! dtransp is in mm d-1 - tile_fluxes(lu)%plant(pft)%dtransp = (1.6 & ! 1.6 - * tile_fluxes(lu)%plant(pft)%gs_accl * tile(lu)%canopy%fapar * climate%dppfd & ! gs - * climate%dvpd) & ! D - * h2o_molmass * (1.0d0 / rho_water) & - * myinterface%params_siml%secs_per_tstep ! convert: mol m-2 s-1 * kg-h2o mol-1 * m3 kg-1 * s day-1 * mm m-1 = mm day-1 - - else - ! Using Phydro gs - tile_fluxes(lu)%plant(pft)%dtransp = out_phydro_inst%e & ! Phydro e is 1.6 gs D - * h2o_molmass * (1.0d0 / rho_water) & - * myinterface%params_siml%secs_per_tstep ! convert: mol m-2 s-1 * kg-h2o mol-1 * m3 kg-1 * s day-1 * mm m-1 = mm day-1 - - ! ~~ This has been moved to waterbal_splash ~~ - ! tile_fluxes(lu)%canopy%dpet_e_soil = out_phydro_inst%le_s_wet & - ! * myinterface%params_siml%secs_per_tstep ! convert: J m-2 s-1 * s day-1 = J m-2 day-1 - - ! print *, "Canopy ET (mm d-1) = ", tile_fluxes(lu)%canopy%daet_canop - ! print *, "Soil LE (J m-2 d-1) = ", climate%dnetrad, tile_fluxes(lu)%canopy%daet_e_soil - end if + ! We plug Pmodel/Phydro-derived gs into T = 1.6gsD + if (.not. use_phydro) then + ! Using P-model gs + ! Note here that stomatal conductance is already normalized by patm (=gs/patm) so E = 1.6 * (gs/patm) * vpd, which is the same as 1.6 gs (vpd/patm) + ! but it is expressed per unit absorbed light, so multiply by PPFD*fapar + ! dtransp is in mm d-1 + tile_fluxes(lu)%plant(pft)%dtransp = (1.6 & ! 1.6 + * tile_fluxes(lu)%plant(pft)%gs_accl * tile(lu)%canopy%fapar * climate%dppfd & ! gs + * climate%dvpd) & ! D + * h2o_molmass * (1.0d0 / rho_water) & + * myinterface%params_siml%secs_per_tstep ! convert: mol m-2 s-1 * kg-h2o mol-1 * m3 kg-1 * s day-1 * mm m-1 = mm day-1 + + ! print*,'tile_fluxes(lu)%plant(pft)%gs_accl ', tile_fluxes(lu)%plant(pft)%gs_accl + ! print*,'tile_fluxes(lu)%plant(pft)%dtransp ', tile_fluxes(lu)%plant(pft)%dtransp + + else + ! Using Phydro gs + tile_fluxes(lu)%plant(pft)%dtransp = out_phydro_inst%e & ! Phydro e is 1.6 gs D + * h2o_molmass * (1.0d0 / rho_water) & + * myinterface%params_siml%secs_per_tstep ! convert: mol m-2 s-1 * kg-h2o mol-1 * m3 kg-1 * s day-1 * mm m-1 = mm day-1 + end if end do pftloop diff --git a/src/plant_pmodel.mod.f90 b/src/plant_pmodel.mod.f90 index 5f2c4115..f387b32a 100644 --- a/src/plant_pmodel.mod.f90 +++ b/src/plant_pmodel.mod.f90 @@ -44,7 +44,7 @@ module md_plant_pmodel integer :: pftno ! canopy - real :: fpc_grid ! fractional projective cover + real :: fpc_grid ! fractional projective cover, sum over all PFTs must add up to 1 (even if there is bare ground, that's treated by fAPAR) real :: lai_ind ! fraction of absorbed photosynthetically active radiation real :: fapar_ind ! fraction of absorbed photosynthetically active radiation real :: acrown ! crown area diff --git a/src/tile_pmodel.mod.f90 b/src/tile_pmodel.mod.f90 index b42701ca..406297d5 100644 --- a/src/tile_pmodel.mod.f90 +++ b/src/tile_pmodel.mod.f90 @@ -540,9 +540,6 @@ subroutine diag_daily( tile, tile_fluxes ) !---------------------------------------------------------------- ! Sum over PFTs to get canopy-level quantities !---------------------------------------------------------------- - ! xxx test - print*,'Reasonable? tile(lu)%plant(:)%fpc_grid ', tile(lu)%plant(:)%fpc_grid - do lu=1,nlu 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) diff --git a/src/waterbal_splash.mod.f90 b/src/waterbal_splash.mod.f90 index 589d0d80..2d3170c7 100644 --- a/src/waterbal_splash.mod.f90 +++ b/src/waterbal_splash.mod.f90 @@ -399,7 +399,7 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g ! Eq. 81, SPLASH 2.0 Documentation tile_fluxes%canopy%daet = (24.0/pi) * (radians(sw * hi) + rx * rw * rv * (dgsin(hn) - dgsin(hi)) + & radians((rx * rw * ru - rx * tile_fluxes%canopy%rnl) * (hn - hi))) ! JAIDEEP FIXME: Technically correct, but for clarity, apply radians to just (hn-hi) ? - tile_fluxes%canopy%daet_e = tile_fluxes%canopy%daet / energy_to_mm ! JAIDEEP FIXME [resolved]: Oops! This is a case where you should use a simple mass-energy conversion, not econ + tile_fluxes%canopy%daet_e = tile_fluxes%canopy%daet / energy_to_mm if (using_pml) then print*,'Warning: simulation parameter use_pml == .true. but not used in combination with SPLASH-AET (using_gs == .false.).' @@ -413,8 +413,8 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g !--------------------------------------------------------- ! Potential soil evaporation !--------------------------------------------------------- - ! potential soil evaporation, not limited by history of P/PET - dpet_soil = (1.0 - fapar) * tile_fluxes%canopy%drn * tile_fluxes%canopy%econ * 1.0 ! 1000 * econ converts energy into mm evaporation + ! potential soil evaporation, not limited by history of P/PET (mm d-1) + dpet_soil = (1.0 - fapar) * tile_fluxes%canopy%drn * tile_fluxes%canopy%econ * 1000.0 ! econ converts energy into mm evaporation !--------------------------------------------------------- ! soil moisture limitation factor @@ -424,8 +424,13 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g ! but a continuous dampening (low pass filter, using dampen_variability()) is applied here instead of a running sum. p_memory = dampen_variability(climate%dprec * secs_per_day, 30.0, p_memory ) pet_memory = dampen_variability(dpet_soil, 30.0, pet_memory ) - p_over_pet_memory = p_memory/(pet_memory + 1e-6) ! corresponds to f in Zhang et al., 2017 Eq. 9, (+ 1e-6) to avoid division by zero - f_soil_aet = max(min(p_over_pet_memory, 1.0), 0.0) + + if (pet_memory > 0.0) then + p_over_pet_memory = p_memory / pet_memory ! corresponds to f in Zhang et al., 2017 Eq. 9, (+ 1e-6) to avoid division by zero + f_soil_aet = max(min(p_over_pet_memory, 1.0), 0.0) + else + f_soil_aet = 1.0 + end if !--------------------------------------------------------- ! Actual soil evaporation (mm d-1 and J d-1) @@ -479,14 +484,9 @@ subroutine calc_et( tile_fluxes, grid, climate, sw, fapar, using_phydro, using_g tile_fluxes%canopy%daet = tile_fluxes%canopy%daet_canop + tile_fluxes%canopy%daet_soil tile_fluxes%canopy%daet_e = tile_fluxes%canopy%daet_e_canop + tile_fluxes%canopy%daet_e_soil - ! print *, "P (mm d-1), PET (mm d-1), P/PET, Avg(P/PET), f_soil_aet = ", (climate%dprec*86400), & - ! tile_fluxes%canopy%dpet_soil, p_over_pet, & - ! p_over_pet_memory, f_soil_aet - ! print *, "Canopy ET (mm d-1, J m-2 d-1) = ", tile_fluxes%canopy%daet_canop, tile_fluxes%canopy%daet_e_canop - ! print *, "Soil ET (mm d-1, J m-2 d-1) = ", tile_fluxes%canopy%daet_soil, tile_fluxes%canopy%daet_e_soil + ! print*,'waterbal: tile_fluxes%canopy%daet_canop, tile_fluxes%canopy%daet_soil ', tile_fluxes%canopy%daet_canop, tile_fluxes%canopy%daet_soil end if - ! xxx debug ! if (splashtest) then diff --git a/vignettes/pmodel_use_newdata.Rmd b/vignettes/pmodel_use_newdata.Rmd index c165c9bf..947367cc 100644 --- a/vignettes/pmodel_use_newdata.Rmd +++ b/vignettes/pmodel_use_newdata.Rmd @@ -330,30 +330,30 @@ output <- rsofun::runread_pmodel_f( p_model_drivers, par = params_modl ) +``` -print( - output$data[[1]] %>% - select(date, gpp, le) %>% - pivot_longer(-date) %>% - mutate(type="pred") %>% - rbind( - p_model_validation$data[[1]] %>% - select(date, gpp, le) %>% - pivot_longer(-date) %>% - mutate(type="obs") - ) %>% - pivot_wider(names_from = type, values_from = value) %>% - group_by(name) %>% - reframe(obs = obs, pred = pred, density=scale(get_density(pred, obs,n=100))) %>% - ggplot(aes(y=obs, x=pred, colour=density)) + - scale_color_viridis_c() + - geom_point(alpha=0.7) + - geom_abline(slope=1, intercept=0, col="red")+ - theme_classic() + - theme(strip.background = element_rect(color = "white", size = 1))+ - facet_wrap(~name, scales = "free", nrow = 1)+ - labs(colour="Density") -) +```{r} +output$data[[1]] %>% + select(date, gpp, le) %>% + pivot_longer(-date) %>% + mutate(type="pred") %>% + rbind( + p_model_validation$data[[1]] %>% + select(date, gpp, le) %>% + pivot_longer(-date) %>% + mutate(type="obs") + ) %>% + pivot_wider(names_from = type, values_from = value) %>% + group_by(name) %>% + reframe(obs = obs, pred = pred, density=scale(get_density(pred, obs,n=100))) %>% + ggplot(aes(y=obs, x=pred, colour=density)) + + scale_color_viridis_c() + + geom_point(alpha=0.7) + + geom_abline(slope=1, intercept=0, col="red")+ + theme_classic() + + theme(strip.background = element_rect(color = "white", size = 1))+ + facet_wrap(~name, scales = "free", nrow = 1)+ + labs(colour="Density") ``` @@ -374,7 +374,7 @@ p_model_drivers <- ) ``` -Run the model and plot outputs +Run the model. ```{r} # run the model for these parameters @@ -382,30 +382,32 @@ output <- rsofun::runread_pmodel_f( p_model_drivers, par = params_modl ) +``` -print( - output$data[[1]] %>% - select(date, gpp, le) %>% - pivot_longer(-date) %>% - mutate(type="pred") %>% - rbind( - p_model_validation$data[[1]] %>% - select(date, gpp, le) %>% - pivot_longer(-date) %>% - mutate(type="obs") - ) %>% - pivot_wider(names_from = type, values_from = value) %>% - group_by(name) %>% - reframe(obs = obs, pred = pred, density=scale(get_density(pred, obs,n=100))) %>% - ggplot(aes(y=obs, x=pred, colour=density)) + - scale_color_viridis_c() + - geom_point(alpha=0.7) + - geom_abline(slope=1, intercept=0, col="red")+ - theme_classic() + - theme(strip.background = element_rect(color = "white", size = 1))+ - facet_wrap(~name, scales = "free", nrow = 1)+ - labs(colour="Density") -) +Plot outputs. + +```{r} +output$data[[1]] %>% + select(date, gpp, le) %>% + pivot_longer(-date) %>% + mutate(type="pred") %>% + rbind( + p_model_validation$data[[1]] %>% + select(date, gpp, le) %>% + pivot_longer(-date) %>% + mutate(type="obs") + ) %>% + pivot_wider(names_from = type, values_from = value) %>% + group_by(name) %>% + reframe(obs = obs, pred = pred, density=scale(get_density(pred, obs,n=100))) %>% + ggplot(aes(y=obs, x=pred, colour=density)) + + scale_color_viridis_c() + + geom_point(alpha=0.7) + + geom_abline(slope=1, intercept=0, col="red")+ + theme_classic() + + theme(strip.background = element_rect(color = "white", size = 1))+ + facet_wrap(~name, scales = "free", nrow = 1)+ + labs(colour="Density") ``` ## Gs-coupled P-hydro run