tested: identical results for gs-coupled ET with diffusion (not PML).

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
 

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

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

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)

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

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,38 +374,40 @@ p_model_drivers <-
       )
 ```
 
-Run the model and plot outputs
+Run the model.
 
 ```{r}
 # run the model for these parameters
 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