diff --git a/ED/build/bin/dependency.mk b/ED/build/bin/dependency.mk index 5783ef47e..d7609f2b8 100644 --- a/ED/build/bin/dependency.mk +++ b/ED/build/bin/dependency.mk @@ -13,7 +13,7 @@ bdf2_solver.o: consts_coms.mod ed_misc_coms.mod ed_state_vars.mod bdf2_solver.o: ed_therm_lib.mod grid_coms.mod rk4_coms.mod soil_coms.mod bdf2_solver.o: therm_lib8.mod canopy_struct_dynamics.o: allometry.mod canopy_air_coms.mod -canopy_struct_dynamics.o: canopy_layer_coms.mod consts_coms.mod +canopy_struct_dynamics.o: canopy_layer_coms.mod consts_coms.mod ed_misc_coms.mod canopy_struct_dynamics.o: ed_state_vars.mod grid_coms.mod met_driver_coms.mod canopy_struct_dynamics.o: pft_coms.mod phenology_coms.mod physiology_coms.mod canopy_struct_dynamics.o: rk4_coms.mod soil_coms.mod therm_lib.mod @@ -32,13 +32,14 @@ events.o: ed_misc_coms.mod ed_state_vars.mod ed_therm_lib.mod events.o: fuse_fiss_utils.mod grid_coms.mod pft_coms.mod therm_lib.mod farq_leuning.o: c34constants.mod consts_coms.mod pft_coms.mod phenology_coms.mod farq_leuning.o: physiology_coms.mod rk4_coms.mod therm_lib8.mod -fire.o: allometry.mod consts_coms.mod disturb_coms.mod ed_misc_coms.mod -fire.o: ed_state_vars.mod grid_coms.mod soil_coms.mod -forestry.o: allometry.mod disturb_coms.mod disturbance_utils.mod ed_max_dims.mod +fire.o: consts_coms.mod disturb_coms.mod ed_misc_coms.mod ed_state_vars.mod +fire.o: grid_coms.mod soil_coms.mod +forestry.o: disturb_coms.mod disturbance_utils.mod ed_max_dims.mod forestry.o: ed_misc_coms.mod ed_state_vars.mod fuse_fiss_utils.mod grid_coms.mod growth_balive.o: allometry.mod consts_coms.mod decomp_coms.mod ed_max_dims.mod growth_balive.o: ed_misc_coms.mod ed_state_vars.mod ed_therm_lib.mod -growth_balive.o: grid_coms.mod mortality.mod pft_coms.mod physiology_coms.mod +growth_balive.o: fuse_fiss_utils.mod grid_coms.mod mortality.mod pft_coms.mod +growth_balive.o: physiology_coms.mod heun_driver.o: canopy_air_coms.mod consts_coms.mod ed_max_dims.mod heun_driver.o: ed_misc_coms.mod ed_state_vars.mod grid_coms.mod heun_driver.o: hydrology_coms.mod met_driver_coms.mod rk4_coms.mod @@ -215,10 +216,9 @@ budget_utils.o: therm_lib.mod dateutils.o: consts_coms.mod ed_filelist.o: ed_max_dims.mod ed_grid.o: consts_coms.mod ed_max_dims.mod ed_node_coms.mod grid_coms.mod -ed_therm_lib.o: allometry.mod canopy_air_coms.mod consts_coms.mod -ed_therm_lib.o: ed_max_dims.mod ed_misc_coms.mod ed_state_vars.mod grid_coms.mod -ed_therm_lib.o: pft_coms.mod rk4_coms.mod soil_coms.mod therm_lib.mod -ed_therm_lib.o: therm_lib8.mod +ed_therm_lib.o: canopy_air_coms.mod consts_coms.mod ed_max_dims.mod +ed_therm_lib.o: ed_misc_coms.mod ed_state_vars.mod grid_coms.mod pft_coms.mod +ed_therm_lib.o: rk4_coms.mod soil_coms.mod therm_lib.mod therm_lib8.mod fatal_error.o: ed_node_coms.mod fuse_fiss_utils.o: allometry.mod canopy_layer_coms.mod decomp_coms.mod fuse_fiss_utils.o: disturb_coms.mod ed_max_dims.mod ed_misc_coms.mod diff --git a/ED/dbgbuild/bin/dependency.mk b/ED/dbgbuild/bin/dependency.mk index 5783ef47e..d7609f2b8 100644 --- a/ED/dbgbuild/bin/dependency.mk +++ b/ED/dbgbuild/bin/dependency.mk @@ -13,7 +13,7 @@ bdf2_solver.o: consts_coms.mod ed_misc_coms.mod ed_state_vars.mod bdf2_solver.o: ed_therm_lib.mod grid_coms.mod rk4_coms.mod soil_coms.mod bdf2_solver.o: therm_lib8.mod canopy_struct_dynamics.o: allometry.mod canopy_air_coms.mod -canopy_struct_dynamics.o: canopy_layer_coms.mod consts_coms.mod +canopy_struct_dynamics.o: canopy_layer_coms.mod consts_coms.mod ed_misc_coms.mod canopy_struct_dynamics.o: ed_state_vars.mod grid_coms.mod met_driver_coms.mod canopy_struct_dynamics.o: pft_coms.mod phenology_coms.mod physiology_coms.mod canopy_struct_dynamics.o: rk4_coms.mod soil_coms.mod therm_lib.mod @@ -32,13 +32,14 @@ events.o: ed_misc_coms.mod ed_state_vars.mod ed_therm_lib.mod events.o: fuse_fiss_utils.mod grid_coms.mod pft_coms.mod therm_lib.mod farq_leuning.o: c34constants.mod consts_coms.mod pft_coms.mod phenology_coms.mod farq_leuning.o: physiology_coms.mod rk4_coms.mod therm_lib8.mod -fire.o: allometry.mod consts_coms.mod disturb_coms.mod ed_misc_coms.mod -fire.o: ed_state_vars.mod grid_coms.mod soil_coms.mod -forestry.o: allometry.mod disturb_coms.mod disturbance_utils.mod ed_max_dims.mod +fire.o: consts_coms.mod disturb_coms.mod ed_misc_coms.mod ed_state_vars.mod +fire.o: grid_coms.mod soil_coms.mod +forestry.o: disturb_coms.mod disturbance_utils.mod ed_max_dims.mod forestry.o: ed_misc_coms.mod ed_state_vars.mod fuse_fiss_utils.mod grid_coms.mod growth_balive.o: allometry.mod consts_coms.mod decomp_coms.mod ed_max_dims.mod growth_balive.o: ed_misc_coms.mod ed_state_vars.mod ed_therm_lib.mod -growth_balive.o: grid_coms.mod mortality.mod pft_coms.mod physiology_coms.mod +growth_balive.o: fuse_fiss_utils.mod grid_coms.mod mortality.mod pft_coms.mod +growth_balive.o: physiology_coms.mod heun_driver.o: canopy_air_coms.mod consts_coms.mod ed_max_dims.mod heun_driver.o: ed_misc_coms.mod ed_state_vars.mod grid_coms.mod heun_driver.o: hydrology_coms.mod met_driver_coms.mod rk4_coms.mod @@ -215,10 +216,9 @@ budget_utils.o: therm_lib.mod dateutils.o: consts_coms.mod ed_filelist.o: ed_max_dims.mod ed_grid.o: consts_coms.mod ed_max_dims.mod ed_node_coms.mod grid_coms.mod -ed_therm_lib.o: allometry.mod canopy_air_coms.mod consts_coms.mod -ed_therm_lib.o: ed_max_dims.mod ed_misc_coms.mod ed_state_vars.mod grid_coms.mod -ed_therm_lib.o: pft_coms.mod rk4_coms.mod soil_coms.mod therm_lib.mod -ed_therm_lib.o: therm_lib8.mod +ed_therm_lib.o: canopy_air_coms.mod consts_coms.mod ed_max_dims.mod +ed_therm_lib.o: ed_misc_coms.mod ed_state_vars.mod grid_coms.mod pft_coms.mod +ed_therm_lib.o: rk4_coms.mod soil_coms.mod therm_lib.mod therm_lib8.mod fatal_error.o: ed_node_coms.mod fuse_fiss_utils.o: allometry.mod canopy_layer_coms.mod decomp_coms.mod fuse_fiss_utils.o: disturb_coms.mod ed_max_dims.mod ed_misc_coms.mod diff --git a/ED/run/ED2IN b/ED/run/ED2IN index 0cfcc9814..40d937f2c 100644 --- a/ED/run/ED2IN +++ b/ED/run/ED2IN @@ -410,6 +410,39 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! ISOILCOL -- LEAF-3 and ED-2 soil colour classes that the model will use when ISOILFLG ! + ! is set to 2. Soil classes are from 1 to 20 (1 = lightest; 20 = darkest). ! + ! The values are the same as CLM-4.0. The table is the albedo for visible ! + ! and near infra-red. ! + !---------------------------------------------------------------------------------------! + ! ! + ! |-----------------------------------------------------------------------| ! + ! | | Dry soil | Saturated | | Dry soil | Saturated | ! + ! | Class |-------------+-------------| Class +-------------+-------------| ! + ! | | VIS | NIR | VIS | NIR | | VIS | NIR | VIS | NIR | ! + ! |-------+------+------+------+------+-------+------+------+------+------| ! + ! | 1 | 0.36 | 0.61 | 0.25 | 0.50 | 11 | 0.24 | 0.37 | 0.13 | 0.26 | ! + ! | 2 | 0.34 | 0.57 | 0.23 | 0.46 | 12 | 0.23 | 0.35 | 0.12 | 0.24 | ! + ! | 3 | 0.32 | 0.53 | 0.21 | 0.42 | 13 | 0.22 | 0.33 | 0.11 | 0.22 | ! + ! | 4 | 0.31 | 0.51 | 0.20 | 0.40 | 14 | 0.20 | 0.31 | 0.10 | 0.20 | ! + ! | 5 | 0.30 | 0.49 | 0.19 | 0.38 | 15 | 0.18 | 0.29 | 0.09 | 0.18 | ! + ! | 6 | 0.29 | 0.48 | 0.18 | 0.36 | 16 | 0.16 | 0.27 | 0.08 | 0.16 | ! + ! | 7 | 0.28 | 0.45 | 0.17 | 0.34 | 17 | 0.14 | 0.25 | 0.07 | 0.14 | ! + ! | 8 | 0.27 | 0.43 | 0.16 | 0.32 | 18 | 0.12 | 0.23 | 0.06 | 0.12 | ! + ! | 9 | 0.26 | 0.41 | 0.15 | 0.30 | 19 | 0.10 | 0.21 | 0.05 | 0.10 | ! + ! | 10 | 0.25 | 0.39 | 0.14 | 0.28 | 20 | 0.08 | 0.16 | 0.04 | 0.08 | ! + ! |-----------------------------------------------------------------------| ! + ! ! + ! Soil type 21 is a special case in which we use the albedo method that used to be ! + ! the default in ED-2.1. ! + !---------------------------------------------------------------------------------------! + NL%ISOILCOL = 21 + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! ! ISOILCOL -- LEAF-3 and ED-2 soil colour classes that the model will use when ISOILFLG ! ! is set to 2. Soil classes are from 1 to 20 (1 = lightest; 20 = darkest). ! @@ -669,6 +702,16 @@ $ED_NL + !---------------------------------------------------------------------------------------! + ! IGRASS -- This controls the dynamics and growth calculation for grasses. A new ! + ! grass scheme is now available where bdead = 0, height is a function of bleaf! + ! and growth happens daily. ALS (3/3/12) ! + ! 0: grasses behave like trees as in ED2.1 (old scheme) ! + ! ! + ! 1: new grass scheme as described above ! + !---------------------------------------------------------------------------------------! + NL%IGRASS = 0 + !---------------------------------------------------------------------------------------! ! IPHEN_SCHEME -- It controls the phenology scheme. Even within each scheme, the ! @@ -1090,7 +1133,7 @@ $ED_NL ! 2. Soil conductivity decreases with depth even for constant soil moisture ! ! , otherwise it is the same as 1. ! !---------------------------------------------------------------------------------------! - NL%IPERCOL = 0 + NL%IPERCOL = 1 !---------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/canopy_struct_dynamics.f90 b/ED/src/dynamics/canopy_struct_dynamics.f90 index f355484b1..d6cd4916e 100644 --- a/ED/src/dynamics/canopy_struct_dynamics.f90 +++ b/ED/src/dynamics/canopy_struct_dynamics.f90 @@ -138,12 +138,16 @@ subroutine canopy_turbulence(cpoly,isi,ipa) , cph2o ! ! intent(in) use soil_coms , only : snow_rough & ! intent(in) , soil_rough ! ! intent(in) + use pft_coms , only : is_grass ! ! intent(in) use therm_lib , only : press2exner & ! function , extheta2temp & ! function , tq2enthalpy ! ! function use allometry , only : h2crownbh & ! function + , size2bl & ! function , dbh2bl ! ! function + use ed_misc_coms , only : igrass ! ! intent(in) use phenology_coms , only : elongf_min ! ! intent(in) + implicit none !----- Arguments --------------------------------------------------------------------! type(polygontype) , target :: cpoly ! Current polygon @@ -771,8 +775,15 @@ subroutine canopy_turbulence(cpoly,isi,ipa) !------ Estimate WAI. ------------------------------------------------------! + if (is_grass(ipft) .and. igrass==1) then + !--use actual leaf mass for grass + waiuse = 0.10 * cpatch%nplant(ico) * cpatch%sla(ico) & + * cpatch%bleaf(ico) + else + !--use dbh for trees waiuse = 0.10 * cpatch%nplant(ico) * cpatch%sla(ico) & * dbh2bl(cpatch%dbh(ico),ipft) + end if !---------------------------------------------------------------------------! @@ -843,7 +854,7 @@ subroutine canopy_turbulence(cpoly,isi,ipa) ! Dry grasses only. Create a pseudo TAI so it won't be a ! ! singularity. ! !------------------------------------------------------------------------! - tai_drygrass = elongf_min * dbh2bl(cpatch%dbh(ico),ipft) + tai_drygrass = elongf_min * size2bl(cpatch%dbh(ico),cpatch%hite(ico),ipft) ladcohort = tai_drygrass / (htopcrown - hbotcrown) !------------------------------------------------------------------------! else @@ -1391,8 +1402,11 @@ subroutine canopy_turbulence8(csite,initp,ipa) , kin_visci8 ! ! intent(in) use soil_coms , only : snow_rough8 & ! intent(in) , soil_rough8 ! ! intent(in) + use pft_coms , only : is_grass ! ! intent(in) use allometry , only : h2crownbh & ! function + , size2bl & ! function , dbh2bl ! ! function + use ed_misc_coms , only : igrass ! ! intent(in) use phenology_coms , only : elongf_min ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -1952,8 +1966,15 @@ subroutine canopy_turbulence8(csite,initp,ipa) !------ Estimate WAI. ------------------------------------------------------! - waiuse = 1.d-1 * initp%nplant(ico) * dble(cpatch%sla(ico)) & - * dble(dbh2bl(cpatch%dbh(ico),ipft)) + if (is_grass(ipft) .and. igrass==1) then + !--use actual leaf mass for grass + waiuse = 1.d-1 * initp%nplant(ico) * dble(cpatch%sla(ico)) & + * dble(cpatch%bleaf(ico)) + else + !--use dbh for trees + waiuse = 1.d-1 * initp%nplant(ico) * dble(cpatch%sla(ico)) & + * dble(dbh2bl(cpatch%dbh(ico),ipft)) + end if !---------------------------------------------------------------------------! @@ -2024,7 +2045,7 @@ subroutine canopy_turbulence8(csite,initp,ipa) ! Dry grasses only. Create a pseudo TAI so it won't be a ! ! singularity. ! !------------------------------------------------------------------------! - tai_drygrass = dble(elongf_min * dbh2bl(cpatch%dbh(ico),ipft)) + tai_drygrass = dble(elongf_min * size2bl(cpatch%dbh(ico),cpatch%hite(ico),ipft)) ladcohort = tai_drygrass / (htopcrown - hbotcrown) !------------------------------------------------------------------------! else diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index e41d7b648..dbfa8f72d 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -815,7 +815,8 @@ subroutine apply_disturbances(cgrid) ,cpatch%pft(ico),cpatch%sla(ico) & ,cpatch%lai(ico),cpatch%wai(ico) & ,cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) + end do csite%area(ipa) = csite%area(ipa) / area_fac csite%age(ipa) = csite%age(ipa) * area_fac @@ -1925,7 +1926,8 @@ subroutine plant_patch(csite,np,mzg,pft,density,ntext_soil,green_leaf_factor ,cpatch%dbh(nc),csite%soil_water(:,np),ntext_soil & ,green_leaf_factor,cpatch%paw_avg(nc),cpatch%elongf(nc) & ,cpatch%phenology_status(nc),cpatch%bleaf(nc) & - ,cpatch%broot(nc),cpatch%bsapwood(nc),cpatch%balive(nc) & + ,cpatch%broot(nc),cpatch%bsapwooda(nc) & + ,cpatch%bsapwoodb(nc),cpatch%balive(nc) & ,cpatch%bstorage(nc)) !------------------------------------------------------------------------------------! @@ -1935,14 +1937,13 @@ subroutine plant_patch(csite,np,mzg,pft,density,ntext_soil,green_leaf_factor call area_indices(cpatch%nplant(nc),cpatch%bleaf(nc),cpatch%bdead(nc) & ,cpatch%balive(nc),cpatch%dbh(nc),cpatch%hite(nc),cpatch%pft(nc) & ,cpatch%sla(nc),cpatch%lai(nc),cpatch%wai(nc),cpatch%crown_area(nc) & - ,cpatch%bsapwood(nc)) + ,cpatch%bsapwooda(nc)) !----- Find the new basal area and above-ground biomass. ----------------------------! - cpatch%basarea(nc) = pio4 * cpatch%dbh(nc) * cpatch%dbh(nc) - cpatch%agb(nc) = ed_biomass(cpatch%bdead(nc),cpatch%balive(nc),cpatch%bleaf(nc) & - ,cpatch%pft(nc),cpatch%hite(nc) ,cpatch%bstorage(nc) & - ,cpatch%bsapwood(nc)) + cpatch%basarea(nc)= pio4 * cpatch%dbh(nc) * cpatch%dbh(nc) + cpatch%agb(nc) = ed_biomass(cpatch%bdead(nc),cpatch%bleaf(nc),cpatch%bsapwooda(nc)& + ,cpatch%pft(nc)) cpatch%leaf_temp(nc) = csite%can_temp(np) cpatch%leaf_temp_pv(nc)=csite%can_temp(np) @@ -1954,7 +1955,7 @@ subroutine plant_patch(csite,np,mzg,pft,density,ntext_soil,green_leaf_factor cpatch%wood_fliq(nc) = 0.0 !----- Because we assigned no water, the internal energy is simply hcap*T. ----------! - call calc_veg_hcap(cpatch%bleaf(nc),cpatch%bdead(nc),cpatch%bsapwood(nc) & + call calc_veg_hcap(cpatch%bleaf(nc),cpatch%bdead(nc),cpatch%bsapwooda(nc) & ,cpatch%nplant(nc),cpatch%pft(nc) & ,cpatch%leaf_hcap(nc),cpatch%wood_hcap(nc)) diff --git a/ED/src/dynamics/events.f90 b/ED/src/dynamics/events.f90 index 9cb8a2511..cfeb49cbb 100644 --- a/ED/src/dynamics/events.f90 +++ b/ED/src/dynamics/events.f90 @@ -287,18 +287,19 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) use ed_state_vars,only: edgrid_g, & edtype,polygontype,sitetype, & patchtype,allocate_patchtype,copy_patchtype,deallocate_patchtype - use pft_coms, only:sla,qsw,q,hgt_min, agf_bs + use pft_coms, only:sla,qsw,q,hgt_min, agf_bs, is_grass use disturbance_utils,only: plant_patch use ed_therm_lib, only: calc_veg_hcap,update_veg_energy_cweh use fuse_fiss_utils, only: terminate_cohorts - use allometry, only : bd2dbh, dbh2h, area_indices, ed_biomass + use allometry, only : bd2dbh, dbh2h, bl2dbh, bl2h, area_indices, ed_biomass use consts_coms, only : pio4 + use ed_misc_coms , only : igrass ! ! intent(in) implicit none real(kind=8),intent(in) :: agb_frac8 real(kind=8),intent(in) :: bgb_frac8 real(kind=8),intent(in) :: fol_frac8 real(kind=8),intent(in) :: stor_frac8 - real :: ialloc,bdead_new,bsw_new,bleaf_new,bfr_new,bstore_new + real :: ialloc,bdead_new,bswa_new,bswb_new,bleaf_new,bfr_new,bstore_new real :: agb_frac,bgb_frac,fol_frac,stor_frac real :: old_leaf_hcap real :: old_wood_hcap @@ -310,9 +311,9 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) type(sitetype),pointer :: csite type(patchtype),pointer :: cpatch - agb_frac = real(agb_frac8) - bgb_frac = real(bgb_frac8) - fol_frac = real(fol_frac8) + agb_frac = real(agb_frac8) + bgb_frac = real(bgb_frac8) + fol_frac = real(fol_frac8) stor_frac = real(stor_frac8) print*,"----------------------------" @@ -347,14 +348,20 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) pft = cpatch%pft(ico) !! calc new pool sizes - ialloc = 1.0 / (1.0 + q(pft) + qsw(pft) * cpatch%hite(ico)) - bdead_new = cpatch%bdead(ico)*& - & (1.0-agb_frac*agf_bs(pft) - bgb_frac*(1.0-agf_bs(pft))) - bsw_new = cpatch%balive(ico) * qsw(pft) * cpatch%hite(ico) *& - & ialloc * (1.0-agb_frac*agf_bs(pft) - bgb_frac*(1.0-agf_bs(pft))) - bstore_new = cpatch%bstorage(ico)*(1.0-stor_frac) - bleaf_new = cpatch%balive(ico) * ialloc *(1.0-fol_frac) - bfr_new = cpatch%balive(ico) * q(pft) * ialloc * (1.0-bgb_frac) + + ialloc = 1.0 / (1.0 + q(pft) + qsw(pft) * cpatch%hite(ico)) + bdead_new = cpatch%bdead(ico) * & + (1.0-agb_frac * agf_bs(pft) - bgb_frac*(1.0-agf_bs(pft))) + bswa_new = cpatch%balive(ico) * qsw(pft) * cpatch%hite(ico) & + * agf_bs(pft) * ialloc * (1.0-agb_frac*agf_bs(pft) & + - bgb_frac*(1.0-agf_bs(pft))) + bswb_new = cpatch%balive(ico) * qsw(pft) * cpatch%hite(ico) & + * (1.-agf_bs(pft)) * ialloc & + * (1.0-agb_frac*agf_bs(pft) - bgb_frac*(1.0-agf_bs(pft))) + + bstore_new = cpatch%bstorage(ico) * (1.0-stor_frac) + bleaf_new = cpatch%balive(ico) * ialloc *(1.0-fol_frac) + bfr_new = cpatch%balive(ico) * q(pft) * ialloc * (1.0-bgb_frac) !! move residual frac to debris/litter pools !! For now assume 100% removal [[needs to be updated]] @@ -363,11 +370,13 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) !! update biomass pools !! [[this needs to be more sophisticated]] - cpatch%balive(ico) = max(0.0,bleaf_new + bfr_new + bsw_new) - cpatch%broot(ico) = max(0.0,bfr_new) - cpatch%bsapwood(ico) = max(0.0,bsw_new) - cpatch%bdead(ico) = max(0.0,bdead_new) - cpatch%bstorage(ico) = max(0.0,bstore_new) + cpatch%balive(ico) = max(0.0,bleaf_new + bfr_new + bswa_new + bswb_new) + cpatch%broot(ico) = max(0.0,bfr_new) + cpatch%bsapwooda(ico) = max(0.0,bswa_new) + cpatch%bsapwoodb(ico) = max(0.0,bswb_new) + cpatch%bdead(ico) = max(0.0,bdead_new) + cpatch%bstorage(ico) = max(0.0,bstore_new) + if(bleaf_new .le. tiny(1.0)) then cpatch%phenology_status(ico) = 2 cpatch%elongf(ico) = 0.0 @@ -379,8 +388,13 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) end if if(cpatch%bdead(ico) .gt. tiny(1.0)) then - cpatch%dbh(ico) = bd2dbh(cpatch%pft(ico), cpatch%bdead(ico)) - cpatch%hite(ico) = dbh2h(cpatch%pft(ico), cpatch%dbh(ico)) + if(is_grass(cpatch%pft(ico)) .and. igrass==1) then + cpatch%dbh(ico) = bl2dbh(cpatch%bleaf(ico),cpatch%pft(ico)) + cpatch%hite(ico) = bl2h (cpatch%bleaf(ico),cpatch%pft(ico)) + else + cpatch%dbh(ico) = bd2dbh(cpatch%pft(ico), cpatch%bdead(ico)) + cpatch%hite(ico) = dbh2h (cpatch%pft(ico), cpatch%dbh(ico)) + end if else cpatch%dbh(ico) = 0.0 cpatch%hite(ico) = 0.0 @@ -391,15 +405,12 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) ,cpatch%balive(ico),cpatch%dbh(ico), cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico), cpatch%lai(ico) & ,cpatch%wai(ico), cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) !----- Update basal area and above-ground biomass. -----------------------! - cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) - + cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !-------------------------------------------------------------------------! ! Here we are leaving all water in the branches and twigs... Do not ! ! worry, if there is any, it will go down through shedding the next ! @@ -408,7 +419,7 @@ subroutine event_harvest(agb_frac8,bgb_frac8,fol_frac8,stor_frac8) old_leaf_hcap = cpatch%leaf_hcap(ico) old_wood_hcap = cpatch%wood_hcap(ico) call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico)) call update_veg_energy_cweh(csite,ipa,ico,old_leaf_hcap,old_wood_hcap) @@ -748,7 +759,8 @@ subroutine event_till(rval8) !! update biomass pools cpatch%balive(ico) = 0.0 cpatch%broot(ico) = 0.0 - cpatch%bsapwood(ico) = 0.0 + cpatch%bsapwooda(ico) = 0.0 + cpatch%bsapwoodb(ico) = 0.0 cpatch%bdead(ico) = 0.0 cpatch%bstorage(ico) = 0.0 cpatch%nplant(ico) = 0.0 @@ -762,9 +774,9 @@ subroutine event_till(rval8) cpatch%wood_water(ico) = 0.0 cpatch%wood_fliq(ico) = 0.0 cpatch%wood_temp(ico) = csite%can_temp(ipa) - call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & - ,cpatch%pft(ico) & + call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & + ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico)) cpatch%leaf_energy(ico) = cpatch%leaf_hcap(ico) * cpatch%leaf_temp(ico) cpatch%wood_energy(ico) = cpatch%wood_hcap(ico) * cpatch%wood_temp(ico) diff --git a/ED/src/dynamics/fire.f90 b/ED/src/dynamics/fire.f90 index 3b6374fdc..f61b8db04 100644 --- a/ED/src/dynamics/fire.f90 +++ b/ED/src/dynamics/fire.f90 @@ -21,7 +21,6 @@ subroutine fire_frequency(cgrid) , fire_smoist_depth & ! intent(in) , k_fire_first & ! intent(in) , fire_parameter ! ! intent(in) - use allometry , only : ed_biomass ! ! function use consts_coms , only : wdns & ! intent(in) , wdnsi & ! intent(in) , day_sec ! ! intent(in) @@ -91,7 +90,9 @@ subroutine fire_frequency(cgrid) ! be defined as the above-ground biomass per unit area. ! !------------------------------------------------------------------------------! cohortloop: do ico = 1,cpatch%ncohorts + fuel = fuel + cpatch%nplant(ico) * cpatch%agb(ico) + end do cohortloop !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 70c217dd8..28fae0685 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -317,7 +317,6 @@ subroutine inventory_mat_forests(cpoly,isi,area_mature_primary,agb_mature_primar , patchtype ! ! structure use disturb_coms , only : plantation_rotation & ! intent(in) , mature_harvest_age ! ! intent(in) - use allometry , only : ed_biomass ! ! function implicit none !----- Arguments -----------------------------------------------------------------------! type(polygontype) , target :: cpoly diff --git a/ED/src/dynamics/growth_balive.f90 b/ED/src/dynamics/growth_balive.f90 index 06d8ff16a..482f30c3f 100644 --- a/ED/src/dynamics/growth_balive.f90 +++ b/ED/src/dynamics/growth_balive.f90 @@ -1,3 +1,4 @@ + !==========================================================================================! !==========================================================================================! module growth_balive @@ -9,7 +10,7 @@ module growth_balive - !=======================================================================================! + !==============================================================1=========================! !=======================================================================================! ! This subroutine will update the alive biomass, and compute the respiration terms ! ! other than leaf respiration. ! @@ -27,6 +28,7 @@ subroutine dbalive_dt(cgrid, tfact) , c2n_storage & ! intent(in) , growth_resp_factor & ! intent(in) , storage_turnover_rate & ! intent(in) + , is_grass & ! intent(in) , phenology ! ! intent(in) use physiology_coms , only : N_plant_lim ! ! intent(in) use grid_coms , only : nzg ! ! intent(in) @@ -35,6 +37,10 @@ subroutine dbalive_dt(cgrid, tfact) use allometry , only : area_indices & ! subroutine , ed_biomass ! ! function use mortality , only : mortality_rates ! ! subroutine + use fuse_fiss_utils , only : sort_cohorts ! ! subroutine + use ed_misc_coms , only : igrass ! ! intent(in) + + implicit none !----- Arguments. -------------------------------------------------------------------! type(edtype) , target :: cgrid @@ -185,9 +191,17 @@ subroutine dbalive_dt(cgrid, tfact) ! Allocate plant carbon balance to balive and bstorage. ! !------------------------------------------------------------------------! balive_in = cpatch%balive(ico) - call alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance & - ,nitrogen_uptake & - ,cpoly%green_leaf_factor(ipft,isi)) + + if(is_grass(ipft).and. igrass==1) then + call alloc_plant_c_balance_grass(csite,ipa,ico,salloc,salloci & + ,carbon_balance,nitrogen_uptake & + ,cpoly%green_leaf_factor(ipft,isi)) + call sort_cohorts(cpatch) + else + call alloc_plant_c_balance(csite,ipa,ico,salloc,salloci & + ,carbon_balance,nitrogen_uptake & + ,cpoly%green_leaf_factor(ipft,isi)) + end if !------------------------------------------------------------------------! @@ -242,13 +256,11 @@ subroutine dbalive_dt(cgrid, tfact) ,cpatch%bdead(ico),cpatch%balive(ico),cpatch%dbh(ico) & ,cpatch%hite(ico) ,cpatch%pft(ico),cpatch%sla(ico) & ,cpatch%lai(ico),cpatch%wai(ico),cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) !----- Update above-ground biomass. -------------------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico),cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !------------------------------------------------------------------------! ! It is likely that biomass has changed, therefore, update ! @@ -257,7 +269,7 @@ subroutine dbalive_dt(cgrid, tfact) old_leaf_hcap = cpatch%leaf_hcap(ico) old_wood_hcap = cpatch%wood_hcap(ico) call calc_veg_hcap(cpatch%bleaf(ico) ,cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico)) call update_veg_energy_cweh(csite,ipa,ico,old_leaf_hcap,old_wood_hcap) @@ -313,8 +325,7 @@ subroutine dbalive_dt_eq_0(cgrid, tfact) use grid_coms , only : nzg ! ! intent(in) use ed_therm_lib , only : calc_veg_hcap & ! function , update_veg_energy_cweh ! ! function - use allometry , only : area_indices & ! subroutine - , ed_biomass ! ! function + use allometry , only : area_indices ! ! subroutine use mortality , only : mortality_rates ! ! subroutine implicit none !----- Arguments. -------------------------------------------------------------------! @@ -382,9 +393,13 @@ subroutine dbalive_dt_eq_0(cgrid, tfact) ,tfact,daily_C_gain,csite%avg_daily_temp(ipa)) !----- Subtract maintenance costs from pools. ---------------------------! - cpatch%balive(ico) = cpatch%balive(ico) - cpatch%bleaf(ico) = cpatch%bleaf(ico) - cpatch%broot(ico) = cpatch%broot(ico) + cpatch%balive(ico) = cpatch%balive(ico) & + - cpatch%leaf_maintenance(ico) & + - cpatch%root_maintenance(ico) + cpatch%bleaf(ico) = cpatch%bleaf(ico) & + - cpatch%leaf_maintenance(ico) + cpatch%broot(ico) = cpatch%broot(ico) & + - cpatch%root_maintenance(ico) cpatch%cb(13,ico) = cpatch%cb(13,ico) & - cpatch%leaf_maintenance(ico) & - cpatch%root_maintenance(ico) @@ -525,6 +540,7 @@ end subroutine dbalive_dt_eq_0 !=======================================================================================! ! This subroutine will transfer some of the stored carbon to balive in order to put ! ! the plant back on allometry. ! + ! ----Is this subroutine ever used??? ALS=== ! !---------------------------------------------------------------------------------------! subroutine transfer_C_from_storage(cpatch,ico,salloc,salloci,nitrogen_uptake & ,N_uptake_pot) @@ -533,6 +549,7 @@ subroutine transfer_C_from_storage(cpatch,ico,salloc,salloci,nitrogen_uptake , c2n_storage & ! intent(in) , c2n_stem & ! intent(in) , q & ! intent(in) + , agf_bs & ! intent(in) , qsw ! ! intent(in) use decomp_coms , only : f_labile ! ! intent(in) use allometry , only : dbh2bl ! ! function @@ -569,7 +586,10 @@ subroutine transfer_C_from_storage(cpatch,ico,salloc,salloci,nitrogen_uptake !----- Compute sapwood and fine root biomass. ---------------------------------------! cpatch%broot(ico) = q(ipft) * cpatch%balive(ico) * salloci - cpatch%bsapwood(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) * salloci + cpatch%bsapwooda(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) * salloci & + * agf_bs(ipft) + cpatch%bsapwoodb(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) * salloci & + *(1.-agf_bs(ipft)) !------------------------------------------------------------------------------------! ! N uptake is required since c2n_leaf < c2n_storage. Units are kgN/plant/day. ! @@ -793,6 +813,8 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance , sla & ! intent(in) , q & ! intent(in) , qsw & ! intent(in) + , agf_bs & ! intent(in) + , r_fract & ! intent(in) , c2n_stem ! ! intent(in) use decomp_coms , only : f_labile ! ! intent(in) use allometry , only : dbh2bl ! ! function @@ -816,16 +838,19 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance real :: old_status real :: delta_bleaf real :: delta_broot - real :: delta_bsapwood + real :: delta_bsapwooda + real :: delta_bsapwoodb real :: available_carbon real :: f_total real :: f_bleaf real :: f_broot - real :: f_bsapwood + real :: f_bsapwooda + real :: f_bsapwoodb real :: f_resp real :: tr_bleaf real :: tr_broot - real :: tr_bsapwood + real :: tr_bsapwooda + real :: tr_bsapwoodb real :: bl logical :: on_allometry logical :: time_to_flush @@ -858,15 +883,18 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance ! plant is drought stress (elongf < 1), we do not allow the plant to get back ! ! to full allometry. ! !------------------------------------------------------------------------------! - bl_max = dbh2bl(cpatch%dbh(ico),ipft) * green_leaf_factor & + !--use dbh for trees + bl_max = dbh2bl(cpatch%dbh(ico),ipft) * green_leaf_factor & * cpatch%elongf(ico) balive_max = dbh2bl(cpatch%dbh(ico),ipft) * salloc * cpatch%elongf(ico) !--- Amount that bleaf, broot, and bsapwood are off allometry -----------------! delta_bleaf = max (0.0, bl_max- cpatch%bleaf(ico)) delta_broot = max (0.0, balive_max * q(ipft) * salloci - cpatch%broot(ico)) - delta_bsapwood = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) * salloci & - - cpatch%bsapwood(ico)) + delta_bsapwooda = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) * salloci & + * agf_bs(ipft) - cpatch%bsapwooda(ico)) + delta_bsapwoodb = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) * salloci & + * (1.-agf_bs(ipft))- cpatch%bsapwoodb(ico)) !------------------------------------------------------------------------------! ! If the available carbon is less than what we need to get back to allometry. ! @@ -876,9 +904,11 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance f_bleaf = delta_bleaf / bl_max f_broot = delta_broot / (balive_max * q(ipft) * salloci ) - f_bsapwood = delta_bsapwood / (balive_max * qsw(ipft) * cpatch%hite(ico) & + f_bsapwooda = delta_bsapwooda / (balive_max * qsw(ipft) * cpatch%hite(ico) & * salloci) - f_total = f_bleaf + f_broot + f_bsapwood + f_bsapwoodb = delta_bsapwoodb / (balive_max * qsw(ipft) * cpatch%hite(ico) & + * salloci) + f_total = f_bleaf + f_broot + f_bsapwooda + f_bsapwoodb !------------------------------------------------------------------------------! ! We only allow transfer from storage to living tissues if there is need ! @@ -887,25 +917,28 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance if (f_total > 0.0) then tr_bleaf = min( delta_bleaf , (f_bleaf/f_total) * available_carbon) tr_broot = min( delta_broot , (f_broot/f_total) * available_carbon) - tr_bsapwood = min( delta_bsapwood, (f_bsapwood/f_total) * available_carbon) + tr_bsapwooda = min( delta_bsapwooda, (f_bsapwooda/f_total) * available_carbon) + tr_bsapwoodb = min( delta_bsapwoodb, (f_bsapwoodb/f_total) * available_carbon) else tr_bleaf = 0. tr_broot = 0. - tr_bsapwood = 0. + tr_bsapwooda = 0. + tr_bsapwoodb = 0. end if !------------------------------------------------------------------------------! cpatch%bleaf(ico) = cpatch%bleaf(ico) + tr_bleaf cpatch%broot(ico) = cpatch%broot(ico) + tr_broot - cpatch%bsapwood(ico) = cpatch%bsapwood(ico) + tr_bsapwood + cpatch%bsapwooda(ico) = cpatch%bsapwooda(ico) + tr_bsapwooda + cpatch%bsapwoodb(ico) = cpatch%bsapwoodb(ico) + tr_bsapwoodb cpatch%balive(ico) = cpatch%bleaf(ico) + cpatch%broot(ico) & - + cpatch%bsapwood(ico) + + cpatch%bsapwooda(ico) + cpatch%bsapwoodb(ico) !----- NPP allocation in diff pools in KgC/m2/day. ----------------------------! cpatch%today_nppleaf(ico) = tr_bleaf * cpatch%nplant(ico) cpatch%today_nppfroot(ico) = tr_broot * cpatch%nplant(ico) - cpatch%today_nppsapwood(ico)= tr_bsapwood * cpatch%nplant(ico) + cpatch%today_nppsapwood(ico)= (tr_bsapwooda + tr_bsapwoodb)* cpatch%nplant(ico) cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) !------------------------------------------------------------------------------! @@ -913,7 +946,8 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance ! -allometry, take that from the carbon balance first, then use some of the ! ! storage if needed be. ! !------------------------------------------------------------------------------! - increment = carbon_balance - tr_bleaf - tr_broot - tr_bsapwood + increment = carbon_balance - tr_bleaf - tr_broot - tr_bsapwooda - tr_bsapwoodb + cpatch%bstorage(ico) = max(0.0, cpatch%bstorage(ico) + increment) !------------------------------------------------------------------------------! @@ -1002,8 +1036,10 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance cpatch%balive(ico) = cpatch%balive(ico) - increment cpatch%bleaf(ico) = cpatch%balive(ico) * salloci * green_leaf_factor cpatch%broot(ico) = cpatch%balive(ico) * q(ipft) * salloci - cpatch%bsapwood(ico) = cpatch%balive(ico) * cpatch%hite(ico) * qsw(ipft) & - * salloci + cpatch%bsapwooda(ico) = cpatch%balive(ico) * cpatch%hite(ico) * qsw(ipft) & + * salloci * agf_bs(ipft) + cpatch%bsapwoodb(ico) = cpatch%balive(ico) * cpatch%hite(ico) * qsw(ipft) & + * salloci * (1.-agf_bs(ipft)) cpatch%phenology_status(ico) = 1 else f_resp = cpatch%today_leaf_resp(ico) & @@ -1021,7 +1057,7 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance end if cpatch%balive(ico) = cpatch%bleaf(ico) + cpatch%broot(ico) & - + cpatch%bsapwood(ico) + + cpatch%bsapwooda(ico) + cpatch%bsapwoodb(ico) end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1049,13 +1085,243 @@ subroutine alloc_plant_c_balance(csite,ipa,ico,salloc,salloci,carbon_balance cpatch%today_nppsapwood(ico) = 0.0 cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) end if - + return end subroutine alloc_plant_c_balance !=======================================================================================! !=======================================================================================! +!=======================================================================================! +! Alternative subroutine for grasses +!=======================================================================================! +subroutine alloc_plant_c_balance_grass(csite,ipa,ico,salloc,salloci,carbon_balance & + ,nitrogen_uptake,green_leaf_factor) + use ed_state_vars , only : sitetype & ! structure + , patchtype ! ! structure + use pft_coms , only : c2n_storage & ! intent(in) + , c2n_leaf & ! intent(in) + , sla & ! intent(in) + , q & ! intent(in) + , qsw & ! intent(in) + , agf_bs & ! intent(in) + , r_fract & ! intent(in) + , hgt_max & ! intent(in) + , c2n_stem ! ! intent(in) + use decomp_coms , only : f_labile ! ! intent(in) + use allometry , only : bl2h & ! function + , h2dbh & ! function + , dbh2h ! ! function + implicit none + !----- Arguments. -------------------------------------------------------------------! + type(sitetype) , target :: csite + integer , intent(in) :: ipa + integer , intent(in) :: ico + real , intent(in) :: salloc + real , intent(in) :: salloci + real , intent(in) :: carbon_balance + real , intent(inout) :: nitrogen_uptake + real , intent(in) :: green_leaf_factor + !----- Local variables. -------------------------------------------------------------! + type(patchtype), pointer :: cpatch + integer :: ipft + real :: bl_max + real :: balive_max + real :: bl_pot + real :: increment + real :: old_status + real :: delta_balive + real :: delta_bleaf + real :: delta_broot + real :: delta_bsapwooda + real :: delta_bsapwoodb + real :: available_carbon + real :: f_total + real :: f_bleaf + real :: f_broot + real :: f_bsapwooda + real :: f_bsapwoodb + real :: f_resp + real :: tr_bleaf + real :: tr_broot + real :: tr_bsapwooda + real :: tr_bsapwoodb + real :: bl + real :: dbh_to_height + logical :: on_allometry + logical :: time_to_flush + !------------------------------------------------------------------------------------! + + cpatch => csite%patch(ipa) + + ipft = cpatch%pft(ico) + + !------------------------------------------------------------------------------------! + ! When plants transit from dormancy to leaf flushing, it is possible that ! + ! carbon_balance is negative, but the sum of carbon_balance and bstorage is ! + ! positive. Under this circumstance, we have to allow plants to grow leaves. ! + !------------------------------------------------------------------------------------! + increment = cpatch%bstorage(ico) + carbon_balance + time_to_flush = (carbon_balance <= 0.0) .and. (increment > 0.0) .and. & + (cpatch%phenology_status(ico) == 1) + + if (carbon_balance > 0.0 .or. time_to_flush) then + if ((cpatch%hite(ico)*(1 + 1e-4)) < hgt_max(ipft)) then ! - could use repro_min_h here instead + !------------------------------------------------------------------------------! + ! The grass is in a vegetative growth phase, put carbon into growth. ! + !------------------------------------------------------------------------------! + !--allow grass to use carbon from that day and from storage to grow + available_carbon = carbon_balance + cpatch%bstorage(ico) + + !--scale maximum growth by elongf (currently grass is "evergreen" so elongf=1) + delta_balive = available_carbon * cpatch%elongf(ico) + increment = available_carbon * (1. - cpatch%elongf(ico)) + + delta_bleaf = delta_balive * salloci * green_leaf_factor + delta_broot = delta_balive * salloci * q (ipft) + delta_bsapwooda = delta_balive * salloci * qsw(ipft) * cpatch%hite(ico) & + * agf_bs(ipft) + delta_bsapwoodb = delta_balive * salloci * qsw(ipft) * cpatch%hite(ico) & + * (1.-agf_bs(ipft)) + + + + cpatch%bleaf(ico) = cpatch%bleaf(ico) + delta_bleaf + cpatch%broot(ico) = cpatch%broot(ico) + delta_broot + cpatch%bsapwooda(ico) = cpatch%bsapwooda(ico) + delta_bsapwooda + cpatch%bsapwoodb(ico) = cpatch%bsapwoodb(ico) + delta_bsapwoodb + + cpatch%balive(ico) = cpatch%bleaf(ico) + cpatch%broot(ico) & + + cpatch%bsapwooda(ico) + cpatch%bsapwoodb(ico) + + !----- NPP allocation in diff pools in KgC/m2/day. ----------------------------! + cpatch%today_nppleaf(ico) = delta_bleaf * cpatch%nplant(ico) + cpatch%today_nppfroot(ico) = delta_broot * cpatch%nplant(ico) + cpatch%today_nppsapwood(ico)= (delta_bsapwooda + delta_bsapwoodb)* cpatch%nplant(ico) + cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) + + !----- update height for grasses to match new leaf mass -----------------------! + cpatch%hite(ico) = min(hgt_max(ipft), bl2h(cpatch%bleaf(ico), ipft)) !limit by maximum height + cpatch%dbh(ico) = h2dbh(cpatch%hite(ico), ipft) !--effective_dbh value for grasses + + + !----- put remaining carbon in the storage pool -------------------------------! + cpatch%bstorage(ico) = max(0.0, cpatch%bstorage(ico) + increment) + !------------------------------------------------------------------------------! + + + if (increment <= 0.0) then + !---------------------------------------------------------------------------! + ! We are using up all of daily C gain and some of bstorage. First ! + ! calculate N demand from using daily C gain. ! + !---------------------------------------------------------------------------! + + nitrogen_uptake = nitrogen_uptake + carbon_balance & + * ( f_labile(ipft) / c2n_leaf(ipft) & + + (1.0 - f_labile(ipft)) / c2n_stem(ipft) ) + + !------------------------------------------------------------------------! + ! Now calculate additional N uptake required from transfer of C from ! + ! storage to balive. ! + !------------------------------------------------------------------------! + nitrogen_uptake = nitrogen_uptake + ( - 1.0 * increment ) & + * ( f_labile(ipft) / c2n_leaf(ipft) & + + (1.0 - f_labile(ipft)) / c2n_stem(ipft) & + - 1.0 / c2n_storage) + + else + !---------------------------------------------------------------------------! + ! N uptake for fraction of daily C gain going to balive. ! + !---------------------------------------------------------------------------! + nitrogen_uptake = nitrogen_uptake + (carbon_balance - increment) & + * ( f_labile(ipft) / c2n_leaf(ipft) & + + (1.0 - f_labile(ipft)) / c2n_stem(ipft)) + !----- N uptake for fraction of daily C gain going to bstorage. ------------! + nitrogen_uptake = nitrogen_uptake + increment / c2n_storage + end if + + else !-- plant is at the maximum height. + !------------------------------------------------------------------------------! + ! Grass is at its maximum height. Put carbon gain into storage. ! + ! For agriculture, put carbon into grain (and storage?) ! + !------------------------------------------------------------------------------! + + !--test here if pft is agriculture, if so put most carbon into grain and ------! + !- maybe a little into storage -- STILL TO BE WRITTEN! ------------------------! + !----- Consider adding allocation to harvest pool here ------------------------! + !---ALS=== Agriculture + increment = carbon_balance ! subtract the part that goes into grain for ag here + + cpatch%bstorage(ico) = cpatch%bstorage(ico) + increment + nitrogen_uptake = nitrogen_uptake + increment / c2n_storage + + + + !----- NPP allocation in diff pools in Kg C/m2/day. ---------------------------! + cpatch%today_nppleaf(ico) = 0.0 + cpatch%today_nppfroot(ico) = 0.0 + cpatch%today_nppsapwood(ico) = 0.0 + cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) + end if !-- end height loop for cb>0 + + + else !-- carbon_balance <0 + !---------------------------------------------------------------------------------! + ! Carbon balance is negative, take it out of storage. ! + !---------------------------------------------------------------------------------! + increment = cpatch%bstorage(ico) + carbon_balance + + + + if (increment <= 0.0) then !-- carbon_balance > bstorage + !----- Take carbon from the Storage pool first --------------------------------! + increment = - increment + cpatch%bstorage(ico) = 0.0 + csite%fsn_in(ipa) = csite%fsn_in(ipa) + cpatch%bstorage(ico) / c2n_storage & + * cpatch%nplant(ico) + + !-- Take the remaining carbon from balive -------------------------------------! + cpatch%balive(ico) = cpatch%balive(ico) - increment + cpatch%bleaf(ico) = cpatch%balive(ico) * salloci * green_leaf_factor + cpatch%broot(ico) = cpatch%balive(ico) * salloci * q(ipft) + cpatch%bsapwooda(ico) = cpatch%balive(ico) * salloci * qsw(ipft) & + * cpatch%hite(ico) * agf_bs(ipft) + cpatch%bsapwoodb(ico) = cpatch%balive(ico) * salloci * qsw(ipft) & + * cpatch%hite(ico) * (1. - agf_bs(ipft)) + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! NOT SURE IF THIS IS CORRECT N ACCOUNTING ! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + csite%fsn_in(ipa) = csite%fsn_in(ipa) + increment & + * ( f_labile(ipft) / c2n_leaf(ipft) & + + (1.0 - f_labile(ipft)) / c2n_stem(ipft) ) & + * cpatch%nplant(ico) + + else !-- carbon_balance < bstorage + !------ Remove the necessary carbon from the storage pool. -------------------! + !------ Dont' forget the nitrogen. --------------------------------------------! + cpatch%bstorage(ico) = cpatch%bstorage(ico) + carbon_balance + csite%fsn_in(ipa) = csite%fsn_in(ipa) - carbon_balance & + * ( f_labile(ipft) / c2n_leaf(ipft) & + + (1.0 - f_labile(ipft)) / c2n_stem(ipft)) & + * cpatch%nplant(ico) + end if + + !---- NPP allocation in diff pools in KgC/m2/day. --------------------------------! + cpatch%today_nppleaf(ico) = 0.0 + cpatch%today_nppfroot(ico) = 0.0 + cpatch%today_nppsapwood(ico) = 0.0 + cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) + end if + + + return +end subroutine alloc_plant_c_balance_grass +!=======================================================================================! +!=======================================================================================! @@ -1074,6 +1340,7 @@ subroutine alloc_plant_c_balance_eq_0(csite,ipa,ico,salloc,salloci,carbon_balanc , sla & ! intent(in) , q & ! intent(in) , qsw & ! intent(in) + , agf_bs & ! intent(in) , c2n_stem ! ! intent(in) use decomp_coms , only : f_labile ! ! intent(in) use allometry , only : dbh2bl ! ! function @@ -1097,16 +1364,19 @@ subroutine alloc_plant_c_balance_eq_0(csite,ipa,ico,salloc,salloci,carbon_balanc real :: old_status real :: delta_bleaf real :: delta_broot - real :: delta_bsapwood + real :: delta_bsapwooda + real :: delta_bsapwoodb real :: available_carbon real :: f_total real :: f_bleaf real :: f_broot - real :: f_bsapwood + real :: f_bsapwooda + real :: f_bsapwoodb real :: f_resp real :: tr_bleaf real :: tr_broot - real :: tr_bsapwood + real :: tr_bsapwooda + real :: tr_bsapwoodb real :: bl logical :: on_allometry logical :: time_to_flush @@ -1146,8 +1416,10 @@ subroutine alloc_plant_c_balance_eq_0(csite,ipa,ico,salloc,salloci,carbon_balanc !--- Amount that bleaf, broot, and bsapwood are off allometry -----------------! delta_bleaf = max (0.0, bl_max- cpatch%bleaf(ico)) delta_broot = max (0.0, balive_max * q(ipft) * salloci - cpatch%broot(ico)) - delta_bsapwood = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) * salloci & - - cpatch%bsapwood(ico)) + delta_bsapwooda = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) & + * salloci * agf_bs(ipft) - cpatch%bsapwooda(ico)) + delta_bsapwoodb = max (0.0, balive_max * qsw(ipft) * cpatch%hite(ico) & + * salloci * (1.-agf_bs(ipft)) - cpatch%bsapwoodb(ico)) !------------------------------------------------------------------------------! ! If the available carbon is less than what we need to get back to allometry. ! @@ -1155,38 +1427,43 @@ subroutine alloc_plant_c_balance_eq_0(csite,ipa,ico,salloc,salloci,carbon_balanc ! extra into bstorage. ! !------------------------------------------------------------------------------! - f_bleaf = delta_bleaf / bl_max - f_broot = delta_broot / (balive_max * q(ipft) * salloci ) - f_bsapwood = delta_bsapwood / (balive_max * qsw(ipft) * cpatch%hite(ico) & - * salloci) - f_total = f_bleaf + f_broot + f_bsapwood + f_bleaf = delta_bleaf / bl_max + f_broot = delta_broot / (balive_max * q(ipft) * salloci ) + f_bsapwooda = delta_bsapwooda / (balive_max * qsw(ipft) * cpatch%hite(ico) & + * salloci * agf_bs(ipft)) + f_bsapwoodb = delta_bsapwoodb / (balive_max * qsw(ipft) * cpatch%hite(ico) & + * salloci * (1.-agf_bs(ipft))) + f_total = f_bleaf + f_broot + f_bsapwooda + f_bsapwoodb !------------------------------------------------------------------------------! ! We only allow transfer from storage to living tissues if there is need ! ! to transfer. ! !------------------------------------------------------------------------------! if (f_total > 0.0) then - tr_bleaf = min( delta_bleaf , (f_bleaf/f_total) * available_carbon) - tr_broot = min( delta_broot , (f_broot/f_total) * available_carbon) - tr_bsapwood = min( delta_bsapwood, (f_bsapwood/f_total) * available_carbon) + tr_bleaf = min(delta_bleaf , (f_bleaf/f_total) * available_carbon) + tr_broot = min(delta_broot , (f_broot/f_total) * available_carbon) + tr_bsapwooda = min(delta_bsapwooda, (f_bsapwooda/f_total) * available_carbon) + tr_bsapwoodb = min(delta_bsapwoodb, (f_bsapwoodb/f_total) * available_carbon) else - tr_bleaf = 0. - tr_broot = 0. - tr_bsapwood = 0. + tr_bleaf = 0. + tr_broot = 0. + tr_bsapwooda = 0. + tr_bsapwoodb = 0. end if !------------------------------------------------------------------------------! - cpatch%bleaf(ico) = cpatch%bleaf(ico) - cpatch%broot(ico) = cpatch%broot(ico) - cpatch%bsapwood(ico) = cpatch%bsapwood(ico) + cpatch%bleaf(ico) = cpatch%bleaf(ico) + tr_bleaf + cpatch%broot(ico) = cpatch%broot(ico) + tr_broot + cpatch%bsapwooda(ico) = cpatch%bsapwooda(ico) + tr_bsapwooda + cpatch%bsapwoodb(ico) = cpatch%bsapwoodb(ico) + tr_bsapwoodb cpatch%balive(ico) = cpatch%bleaf(ico) + cpatch%broot(ico) & - + cpatch%bsapwood(ico) + + cpatch%bsapwooda(ico) + cpatch%bsapwoodb(ico) !----- NPP allocation in diff pools in KgC/m2/day. ----------------------------! - cpatch%today_nppleaf(ico) = tr_bleaf * cpatch%nplant(ico) - cpatch%today_nppfroot(ico) = tr_broot * cpatch%nplant(ico) - cpatch%today_nppsapwood(ico)= tr_bsapwood * cpatch%nplant(ico) + cpatch%today_nppleaf(ico) = tr_bleaf * cpatch%nplant(ico) + cpatch%today_nppfroot(ico) = tr_broot * cpatch%nplant(ico) + cpatch%today_nppsapwood(ico)= (tr_bsapwooda + tr_bsapwoodb) * cpatch%nplant(ico) cpatch%today_nppdaily(ico) = carbon_balance * cpatch%nplant(ico) !------------------------------------------------------------------------------! @@ -1194,7 +1471,7 @@ subroutine alloc_plant_c_balance_eq_0(csite,ipa,ico,salloc,salloci,carbon_balanc ! -allometry, take that from the carbon balance first, then use some of the ! ! storage if needed be. ! !------------------------------------------------------------------------------! - increment = carbon_balance - tr_bleaf - tr_broot - tr_bsapwood + increment = carbon_balance - tr_bleaf - tr_broot - tr_bsapwooda - tr_bsapwoodb cpatch%bstorage(ico) = max(0.0, cpatch%bstorage(ico) + increment) !------------------------------------------------------------------------------! @@ -1294,7 +1571,6 @@ end subroutine alloc_plant_c_balance_eq_0 - !=======================================================================================! !=======================================================================================! subroutine potential_N_uptake(cpatch,ico,salloc,salloci,balive_in,carbon_balance_pot & @@ -1330,7 +1606,7 @@ subroutine potential_N_uptake(cpatch,ico,salloc,salloci,balive_in,carbon_balance N_uptake_pot = N_uptake_pot + carbon_balance_pot / c2n_storage elseif (cpatch%phenology_status(ico) == 1) then - + ! this calculation of bl_max is wrong for grass, but they should not have phenology_status=1 yet bl_max = dbh2bl(cpatch%dbh(ico),ipft) * green_leaf_factor * cpatch%elongf(ico) bl_pot = cpatch%bleaf(ico) + carbon_balance_pot diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index f91e520f0..3c037f229 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -57,6 +57,7 @@ subroutine mortality_rates(cpatch,ipa,ico,avg_daily_temp, patch_age) expmort = max( lnexp_min, min( lnexp_max & , mort2(ipft) * cpatch%cbr_bar(ico))) cpatch%mort_rate(2,ico) = mort1(ipft) / (1. + exp(expmort)) + !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/phenology_aux.f90 b/ED/src/dynamics/phenology_aux.f90 index 3d3c43229..8b72de86c 100644 --- a/ED/src/dynamics/phenology_aux.f90 +++ b/ED/src/dynamics/phenology_aux.f90 @@ -401,8 +401,8 @@ subroutine first_phenology(cgrid) ,cpatch%paw_avg(ico),cpatch%elongf(ico) & ,cpatch%phenology_status(ico) & ,cpatch%bleaf(ico),cpatch%broot(ico) & - ,cpatch%bsapwood(ico),cpatch%balive(ico) & - ,cpatch%bstorage(ico)) + ,cpatch%bsapwooda(ico),cpatch%bsapwoodb(ico)& + ,cpatch%balive(ico),cpatch%bstorage(ico)) !------------------------------------------------------------------------! @@ -411,13 +411,13 @@ subroutine first_phenology(cgrid) ,cpatch%balive(ico),cpatch%dbh(ico),cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & ,cpatch%wai(ico),cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) !------------------------------------------------------------------------! !----- Find heat capacity and vegetation internal energy. ---------------! call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico) ) cpatch%leaf_energy(ico) = cpatch%leaf_hcap(ico) * cpatch%leaf_temp(ico) @@ -456,7 +456,7 @@ end subroutine first_phenology !---------------------------------------------------------------------------------------! subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_soil & ,green_leaf_factor,paw_avg,elongf,phenology_status & - ,bleaf,broot,bsapwood,balive,bstorage) + ,bleaf,broot,bsapwooda,bsapwoodb,balive,bstorage) use soil_coms , only : soil & ! intent(in), look-up table , slz & ! intent(in) , slzt & ! intent(in) @@ -464,10 +464,11 @@ subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_s use phenology_coms, only : spot_phen & ! intent(in) , elongf_min ! ! intent(in) use pft_coms , only : phenology & ! intent(in) + , agf_bs & ! intent(in) , q & ! intent(in) , qsw ! ! intent(in) use ed_max_dims , only : n_pft ! ! intent(in) - use allometry , only : dbh2bl & ! function + use allometry , only : size2bl & ! function , h2crownbh ! ! function implicit none !----- Arguments --------------------------------------------------------------------! @@ -484,7 +485,8 @@ subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_s integer , intent(out) :: phenology_status ! phenology Flag real , intent(out) :: bleaf ! Leaf biomass real , intent(out) :: broot ! Root biomass - real , intent(out) :: bsapwood ! Sapwood biomass + real , intent(out) :: bsapwooda ! Sapwood biomass above ground + real , intent(out) :: bsapwoodb ! Sapwood biomass below ground real , intent(out) :: balive ! Living tissue biomass real , intent(out) :: bstorage ! Storage biomass !----- Local variables --------------------------------------------------------------! @@ -499,10 +501,6 @@ subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_s real :: psi_crit ! Critical point potential real :: mcheight ! Mid-crown height !------------------------------------------------------------------------------------! - - - - !------------------------------------------------------------------------------------! ! Here we decide how to compute the mean available water fraction. ! !------------------------------------------------------------------------------------! @@ -520,10 +518,10 @@ subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_s / (soil(nsoil)%soilwp / soil(nsoil)%slmsts) ** soil(nsoil)%slbs psi_crit = soil(nsoil)%slpots & / (soil(nsoil)%soilld / soil(nsoil)%slmsts) ** soil(nsoil)%slbs - paw_avg = paw_avg + max(0.0, (psi_layer - psi_wilt)) * dslz(k) & / (psi_crit - psi_wilt) end do + paw_avg = paw_avg / abs(slz(kroot)) else !----- Use soil moisture (mass) to determine phenology. --------------------------! @@ -569,12 +567,13 @@ subroutine pheninit_balive_bstorage(mzg,ipft,kroot,height,dbh,soil_water,ntext_s !----- Compute the biomass of living tissues. ---------------------------------------! salloc = 1.0 + q(ipft) + qsw(ipft) * height salloci = 1.0 / salloc - bleaf_max = dbh2bl(dbh,ipft) + bleaf_max = size2bl(dbh,height,ipft) balive_max = bleaf_max * salloc bleaf = bleaf_max * elongf broot = balive_max * q(ipft) * salloci - bsapwood = balive_max * qsw(ipft) * height * salloci - balive = bleaf + broot + bsapwood + bsapwooda = balive_max * qsw(ipft) * height * salloci * agf_bs(ipft) + bsapwoodb = balive_max * qsw(ipft) * height * salloci * (1.-agf_bs(ipft)) + balive = bleaf + broot + bsapwooda + bsapwoodb !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/phenology_driv.f90 b/ED/src/dynamics/phenology_driv.f90 index f6d26673d..879ddb307 100644 --- a/ED/src/dynamics/phenology_driv.f90 +++ b/ED/src/dynamics/phenology_driv.f90 @@ -183,6 +183,7 @@ subroutine update_phenology(doy, cpoly, isi, lat) use ed_misc_coms , only : current_time ! ! intent(in) use allometry , only : area_indices & ! subroutine , ed_biomass & ! function + , size2bl & ! function , dbh2bl ! ! function use phenology_aux , only : daylength ! ! function implicit none @@ -429,9 +430,9 @@ subroutine update_phenology(doy, cpoly, isi, lat) ! 2. The plant has no leaves, but the soil has started to come back to more ! ! moist conditions. Given this situation, leaves can start growing again. ! !------------------------------------------------------------------------------! + cpatch%elongf(ico) = max(0.0, min (1.0, cpatch%paw_avg(ico))) - bl_max = cpatch%elongf(ico) * dbh2bl(cpatch%dbh(ico),ipft) - + bl_max = cpatch%elongf(ico) * size2bl(cpatch%dbh(ico),cpatch%hite(ico),ipft) !----- In case it is too dry, drop all the leaves... --------------------------! if (cpatch%elongf(ico) < elongf_min) then @@ -494,13 +495,11 @@ subroutine update_phenology(doy, cpoly, isi, lat) call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico),cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & - ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Update above-ground biomass. ----------------------------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !---------------------------------------------------------------------------------! ! The leaf biomass of the cohort has changed, update the vegetation energy - ! @@ -508,7 +507,7 @@ subroutine update_phenology(doy, cpoly, isi, lat) !---------------------------------------------------------------------------------! old_leaf_hcap = cpatch%leaf_hcap(ico) old_wood_hcap = cpatch%wood_hcap(ico) - call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwood(ico) & + call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwooda(ico) & ,cpatch%nplant(ico),cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico) ) call update_veg_energy_cweh(csite,ipa,ico,old_leaf_hcap,old_wood_hcap) @@ -828,18 +827,15 @@ subroutine update_phenology_eq_0(doy, cpoly, isi, lat) end select !---------------------------------------------------------------------------------! - !----- Update LAI, WAI, and CAI accordingly. -------------------------------------! call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico),cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & - ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Update above-ground biomass. ----------------------------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !---------------------------------------------------------------------------------! ! The leaf biomass of the cohort has changed, update the vegetation energy - ! @@ -847,7 +843,7 @@ subroutine update_phenology_eq_0(doy, cpoly, isi, lat) !---------------------------------------------------------------------------------! old_leaf_hcap = cpatch%leaf_hcap(ico) old_wood_hcap = cpatch%wood_hcap(ico) - call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwood(ico) & + call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwooda(ico) & ,cpatch%nplant(ico),cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico) ) call update_veg_energy_cweh(csite,ipa,ico,old_leaf_hcap,old_wood_hcap) diff --git a/ED/src/dynamics/photosyn_driv.f90 b/ED/src/dynamics/photosyn_driv.f90 index ffc19eb11..87b4da455 100644 --- a/ED/src/dynamics/photosyn_driv.f90 +++ b/ED/src/dynamics/photosyn_driv.f90 @@ -31,6 +31,7 @@ subroutine canopy_photosynthesis(csite,cmet,mzg,ipa,lsl,ntext_soil use phenology_coms , only : llspan_inf ! ! intent(in) use farq_leuning , only : lphysiol_full ! ! sub-routine use allometry , only : h2crownbh ! ! function + implicit none !----- Arguments -----------------------------------------------------------------------! type(sitetype) , target :: csite ! Current site diff --git a/ED/src/dynamics/reproduction.f90 b/ED/src/dynamics/reproduction.f90 index b0215d496..f7b960777 100644 --- a/ED/src/dynamics/reproduction.f90 +++ b/ED/src/dynamics/reproduction.f90 @@ -44,6 +44,8 @@ subroutine reproduction(cgrid, month) use allometry , only : dbh2bd & ! function , dbh2bl & ! function , h2dbh & ! function + , size2bl & ! function + , dbh2h & ! function , ed_biomass & ! function , area_indices & ! subroutine , dbh2krdepth ! ! function @@ -154,7 +156,6 @@ subroutine reproduction(cgrid, month) !---- This time we loop over PFTs, not cohorts. ----------------------------! pftloop: do ipft = 1, n_pft - !------------------------------------------------------------------------! ! Check to make sure we are including the PFT and that it is not too ! ! cold. ! @@ -168,7 +169,6 @@ subroutine reproduction(cgrid, month) ! this PFT to be in an agriculture patch. ! !---------------------------------------------------------------------! if(csite%dist_type(ipa) /= 1 .or. include_pft_ag(ipft)) then - !------------------------------------------------------------------! ! We assign the recruit in the temporary recruitment structure. ! !------------------------------------------------------------------! @@ -177,6 +177,8 @@ subroutine reproduction(cgrid, month) rectest%wood_temp = csite%can_temp(ipa) rectest%leaf_temp_pv=csite%can_temp(ipa) rectest%wood_temp_pv=csite%can_temp(ipa) + + !- recruits start at minimum height and dbh and bleaf are calculated from that rectest%hite = hgt_min(ipft) rectest%dbh = h2dbh(rectest%hite, ipft) rectest%krdepth = dbh2krdepth(rectest%hite,rectest%dbh & @@ -191,8 +193,8 @@ subroutine reproduction(cgrid, month) ,rectest%paw_avg,rectest%elongf & ,rectest%phenology_status & ,rectest%bleaf,rectest%broot & - ,rectest%bsapwood,rectest%balive & - ,rectest%bstorage) + ,rectest%bsapwooda,rectest%bsapwoodb & + ,rectest%balive,rectest%bstorage) rectest%nplant = csite%repro(ipft,ipa) & / ( rectest%balive + rectest%bdead & @@ -227,9 +229,10 @@ subroutine reproduction(cgrid, month) else !------------------------------------------------------------------! ! If we have reached this branch, we are in an agricultural ! - ! patch. Send the seed litter to the soil pools for ! - ! decomposition. ! + ! patch. Send the seed litter to the soil pools for decomposition.! !------------------------------------------------------------------! + !---ALS=== dont send all seeds to litter! Keep it for harvesting? ! + csite%fast_soil_N(ipa) = csite%fast_soil_N(ipa) & + csite%repro(ipft,ipa) / c2n_recruit(ipft) csite%fast_soil_C(ipa) = csite%fast_soil_C(ipa) & @@ -300,7 +303,8 @@ subroutine reproduction(cgrid, month) cpatch%phenology_status(ico) = recruit(inew)%phenology_status cpatch%bleaf (ico) = recruit(inew)%bleaf cpatch%broot (ico) = recruit(inew)%broot - cpatch%bsapwood (ico) = recruit(inew)%bsapwood + cpatch%bsapwooda (ico) = recruit(inew)%bsapwooda + cpatch%bsapwoodb (ico) = recruit(inew)%bsapwoodb cpatch%balive (ico) = recruit(inew)%balive cpatch%bstorage (ico) = recruit(inew)%bstorage cpatch%leaf_temp (ico) = recruit(inew)%leaf_temp @@ -310,7 +314,6 @@ subroutine reproduction(cgrid, month) !---------------------------------------------------------------------! - !----- Initialise the next variables with zeroes... ------------------! cpatch%leaf_water(ico) = 0.0 cpatch%leaf_fliq (ico) = 0.0 @@ -323,11 +326,9 @@ subroutine reproduction(cgrid, month) ! Computing initial AGB and Basal Area. Their derivatives will be ! ! zero. ! !---------------------------------------------------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) & - ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) + cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) cpatch%dagb_dt(ico) = 0.0 cpatch%dba_dt(ico) = 0.0 @@ -346,10 +347,10 @@ subroutine reproduction(cgrid, month) ,cpatch%bdead(ico),cpatch%balive(ico) & ,cpatch%dbh(ico),cpatch%hite(ico),cpatch%pft(ico) & ,cpatch%sla(ico),cpatch%lai(ico),cpatch%wai(ico) & - ,cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Find heat capacity and vegetation internal energy. ------------! call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico)) @@ -490,10 +491,10 @@ subroutine reproduction(cgrid, month) ,cpatch%bdead(ico),cpatch%balive(ico),cpatch%dbh(ico)& ,cpatch%hite(ico), cpatch%pft(ico),cpatch%sla(ico) & ,cpatch%lai(ico),cpatch%wai(ico) & - ,cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Find heat capacity and vegetation internal energy. ---------! call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%nplant(ico) & + ,cpatch%bsapwooda(ico),cpatch%nplant(ico) & ,cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico)) @@ -587,7 +588,6 @@ subroutine reproduction_eq_0(cgrid, month) use allometry , only : dbh2bd & ! function , dbh2bl & ! function , h2dbh & ! function - , ed_biomass & ! function , area_indices ! ! subroutine use grid_coms , only : nzg ! ! intent(in) implicit none @@ -819,7 +819,7 @@ subroutine seed_dispersal(cpoly,late_spring) nseed_maygo = 0. end if !---------------------------------------------------------------------------! - + !---------------------------------------------------------------------------! ! Spread the seedlings across all patches in this polygon. ! !---------------------------------------------------------------------------! @@ -838,8 +838,6 @@ subroutine seed_dispersal(cpoly,late_spring) * cpoly%area(donsi) !---------------------------------------------------------------------! - - !---------------------------------------------------------------------! ! Include the local dispersal if this is the donor patch. ! !---------------------------------------------------------------------! @@ -848,7 +846,7 @@ subroutine seed_dispersal(cpoly,late_spring) + nseed_stays end if !---------------------------------------------------------------------! - + end do recpaloop2 !------------------------------------------------------------------------! end do recsiloop2 diff --git a/ED/src/dynamics/rk4_misc.f90 b/ED/src/dynamics/rk4_misc.f90 index b99fe3680..c4545f047 100644 --- a/ED/src/dynamics/rk4_misc.f90 +++ b/ED/src/dynamics/rk4_misc.f90 @@ -3433,13 +3433,14 @@ subroutine print_csiteipa(csite, ipa) write (unit=*,fmt='(80a)') ('-',k=1,80) write (unit=*,fmt='(2(a7,1x),8(a12,1x))') & ' PFT','KRDEPTH',' NPLANT',' WAI',' DBH',' BDEAD' & - ,' BSAPWOOD',' WOOD_ENERGY',' WOOD_TEMP',' WOOD_WATER' + ,' BSAPWOODA',' BSAPWOODB',' WOOD_ENERGY',' WOOD_TEMP' & + ,' WOOD_WATER' do ico = 1,cpatch%ncohorts if (cpatch%wood_resolvable(ico)) then write(unit=*,fmt='(2(i7,1x),8(es12.4,1x))') cpatch%pft(ico), cpatch%krdepth(ico) & ,cpatch%nplant(ico),cpatch%wai(ico),cpatch%dbh(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico),cpatch%wood_energy(ico),cpatch%wood_temp(ico) & - ,cpatch%wood_water(ico) + ,cpatch%bsapwooda(ico),cpatch%bsapwoodb(ico),cpatch%wood_energy(ico) & + ,cpatch%wood_temp(ico),cpatch%wood_water(ico) end if end do write (unit=*,fmt='(a)' ) ' ' @@ -3676,12 +3677,12 @@ subroutine print_rk4patch(y,csite,ipa) write (unit=*,fmt='(80a)') ('-',k=1,80) write (unit=*,fmt='(2(a7,1x),5(a12,1x))') & ' PFT','KRDEPTH',' NPLANT',' HEIGHT',' DBH',' BDEAD' & - ,' BSAPWOOD' + ,' BSAPWOODA',' BSAPWOODB' do ico = 1,cpatch%ncohorts if (cpatch%wood_resolvable(ico)) then write(unit=*,fmt='(2(i7,1x),5(es12.4,1x))') cpatch%pft(ico), cpatch%krdepth(ico) & ,cpatch%nplant(ico),cpatch%hite(ico),cpatch%dbh(ico),cpatch%bdead(ico) & - ,cpatch%bsapwood(ico) + ,cpatch%bsapwooda(ico),cpatch%bsapwoodb(ico) end if end do write (unit=*,fmt='(80a)') ('-',k=1,80) diff --git a/ED/src/dynamics/structural_growth.f90 b/ED/src/dynamics/structural_growth.f90 index 1dd84f809..6302f3649 100644 --- a/ED/src/dynamics/structural_growth.f90 +++ b/ED/src/dynamics/structural_growth.f90 @@ -19,6 +19,7 @@ subroutine structural_growth(cgrid, month) , c2n_stem & ! intent(in) , l2n_stem & ! intent(in) , negligible_nplant & ! intent(in) + , is_grass & ! intent(in) , agf_bs ! ! intent(in) use decomp_coms , only : f_labile ! ! intent(in) use ed_max_dims , only : n_pft & ! intent(in) @@ -26,6 +27,8 @@ subroutine structural_growth(cgrid, month) use ed_misc_coms , only : ibigleaf ! ! intent(in) use ed_therm_lib , only : calc_veg_hcap & ! function , update_veg_energy_cweh ! ! function + use ed_misc_coms , only : igrass ! ! intent(in) + implicit none !----- Arguments -----------------------------------------------------------------------! type(edtype) , target :: cgrid @@ -46,6 +49,7 @@ subroutine structural_growth(cgrid, month) real :: salloci real :: balive_in real :: bdead_in + real :: bleaf_in real :: hite_in real :: dbh_in real :: nplant_in @@ -91,6 +95,7 @@ subroutine structural_growth(cgrid, month) !----- Remember inputs in order to calculate increments later on. ----------! balive_in = cpatch%balive(ico) bdead_in = cpatch%bdead(ico) + bleaf_in = cpatch%bleaf(ico) hite_in = cpatch%hite(ico) dbh_in = cpatch%dbh(ico) nplant_in = cpatch%nplant(ico) @@ -126,6 +131,13 @@ subroutine structural_growth(cgrid, month) !----- Grow plants; bdead gets fraction f_bdead of bstorage. ---------------! cpatch%bdead(ico) = cpatch%bdead(ico) + f_bdead * cpatch%bstorage(ico) + if (ibigleaf == 0 ) then + !------ NPP allocation to wood and course roots in KgC /m2 --------------! + cpatch%today_NPPwood(ico) = agf_bs(ipft)*f_bdead*cpatch%bstorage(ico) & + * cpatch%nplant(ico) + cpatch%today_NPPcroot(ico) = (1. - agf_bs(ipft)) * f_bdead & + * cpatch%bstorage(ico) * cpatch%nplant(ico) + end if if (ibigleaf == 0 ) then !------ NPP allocation to wood and course roots in KgC /m2 --------------! @@ -152,8 +164,13 @@ subroutine structural_growth(cgrid, month) cpatch%today_NPPseeds(ico) = f_bseeds * cpatch%bstorage(ico) & * cpatch%nplant(ico) + !---------------------------------------------------------------------------! + ! ALS. If agriculture: set seedling_mortality very low or zero ! + ! to keep all of the seeds for harvest later in the season ! + !---------------------------------------------------------------------------! seed_litter = cpatch%bseeds(ico) * cpatch%nplant(ico) & * seedling_mortality(ipft) + !---------------------------------------------------------------------------! ! Rebalance the plant nitrogen uptake considering the actual alloc- ! @@ -196,7 +213,7 @@ subroutine structural_growth(cgrid, month) !---------------------------------------------------------------------------! old_leaf_hcap = cpatch%leaf_hcap(ico) old_wood_hcap = cpatch%wood_hcap(ico) - call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwood(ico) & + call calc_veg_hcap(cpatch%bleaf(ico),cpatch%bdead(ico),cpatch%bsapwooda(ico)& ,cpatch%nplant(ico),cpatch%pft(ico) & ,cpatch%leaf_hcap(ico),cpatch%wood_hcap(ico) ) call update_veg_energy_cweh(csite,ipa,ico,old_leaf_hcap,old_wood_hcap) @@ -222,15 +239,23 @@ subroutine structural_growth(cgrid, month) !----- Compute the relative carbon balance. --------------------------------! cb_act = 0.0 cb_max = 0.0 - do imonth = 1,12 - cb_act = cb_act + cpatch%cb(imonth,ico) - cb_max = cb_max + cpatch%cb_max(imonth,ico) - end do + + if (is_grass(ipft).and. igrass==1) then !!Grass loop, use past month's carbon balance only + cb_act = cpatch%cb(update_month,ico) + cb_max = cpatch%cb_max(update_month,ico) + else !!Tree loop, use annual average carbon balance + do imonth = 1,12 + cb_act = cb_act + cpatch%cb(imonth,ico) + cb_max = cb_max + cpatch%cb_max(imonth,ico) + end do + end if + if(cb_max > 0.0)then cpatch%cbr_bar(ico) = cb_act / cb_max else cpatch%cbr_bar(ico) = 0.0 end if + !---------------------------------------------------------------------------! @@ -503,8 +528,10 @@ subroutine plant_structural_allocation(ipft,hite,dbh,lat,phen_status,f_bseeds,f_ , r_fract & ! intent(in) , st_fract & ! intent(in) , dbh_crit & ! intent(in) + , hgt_max & ! intent(in) , is_grass ! ! intent(in) - use ed_misc_coms , only : current_time ! ! intent(in) + use ed_misc_coms , only : current_time & ! intent(in) + , igrass ! ! intent(in) use ed_misc_coms , only : ibigleaf ! ! intent(in) implicit none !----- Arguments -----------------------------------------------------------------------! @@ -549,7 +576,6 @@ subroutine plant_structural_allocation(ipft,hite,dbh,lat,phen_status,f_bseeds,f_ !---------------------------------------------------------------------------------------! - select case (ibigleaf) case (0) !------------------------------------------------------------------------------------! @@ -560,21 +586,39 @@ subroutine plant_structural_allocation(ipft,hite,dbh,lat,phen_status,f_bseeds,f_ ! dropping leaves or off allometry. ! !------------------------------------------------------------------------------------! if ((phenology(ipft) /= 2 .or. late_spring) .and. phen_status == 0) then - if (is_grass(ipft) .and. dbh >= dbh_crit(ipft)) then - !---------------------------------------------------------------------------------! - ! Grasses have reached the maximum height, stop growing in size and send ! - ! everything to reproduction. ! - !---------------------------------------------------------------------------------! - f_bseeds = 1.0 - st_fract(ipft) - elseif (hite <= repro_min_h(ipft)) then - !----- The plant is too short, invest as much as it can in growth. ---------------! + !------------------------------------------------------------------------------! + !---ALS=== This is where allocation to seeds is occuring. It will need to be ! + ! modified but I'm leaving it for later --- GRASSES! Want to add a functional ! + ! form to constrain this throughout the season - also consider moving this to ! + ! growth_balive since it isn't actually structural growth ! + !------------------------------------------------------------------------------! + if (is_grass(ipft).and. igrass==1) then !-----------------Grass loop--------------! + if ((hite * (1 + 1.0e-4)) >= hgt_max(ipft)) then + !--------------------------------------------------------------------------! + ! Grasses have reached the maximum height, stop growing in size and send ! + ! everything to reproduction. ! + !--------------------------------------------------------------------------! + f_bseeds = 1.0 - st_fract(ipft) + f_bdead = 0.0 + elseif ((hite * (1 + epsilon(1.))) <= repro_min_h(ipft)) then + !----- The plant is too short, invest as much as it can in growth. --------! + f_bseeds = 0.0 + f_bdead = 0.0 + else ! repro_min_h < hite< hgt_max + !----- Plant is with a certain height, use prescribed reproduction rate. --! + f_bseeds = r_fract(ipft) + f_bdead = 0.0 + end if + elseif (hite <= repro_min_h(ipft)) then !----------------- Tree Loop -------------! + !----- The plant is too short, invest as much as it can in growth. ------------! f_bseeds = 0.0 + f_bdead = 1.0 - st_fract(ipft) - f_bseeds else - !----- Plant is with a certain height, use prescribed reproduction rate. ---------! + !----- Plant is with a certain height, use prescribed reproduction rate. ------! f_bseeds = r_fract(ipft) - end if - f_bdead = 1.0 - st_fract(ipft) - f_bseeds - else + f_bdead = 1.0 - st_fract(ipft) - f_bseeds + end if !!end tree loop + else !-- Plant should not allocate carbon to seeds or grow new biomass -------------! f_bdead = 0.0 f_bseeds = 0.0 end if @@ -611,7 +655,7 @@ subroutine plant_structural_allocation(ipft,hite,dbh,lat,phen_status,f_bseeds,f_ close (unit=66,status='keep') end if !---------------------------------------------------------------------------------------! - + return end subroutine plant_structural_allocation !==========================================================================================! @@ -632,14 +676,19 @@ subroutine update_derived_cohort_props(cpatch,ico,green_leaf_factor,lsl) use ed_state_vars , only : patchtype ! ! structure use pft_coms , only : phenology & ! intent(in) , q & ! intent(in) - , qsw ! ! intent(in) + , qsw & ! intent(in) + , is_grass ! ! intent(in) use allometry , only : bd2dbh & ! function , dbh2h & ! function , dbh2bl & ! function , dbh2krdepth & ! function + , bl2dbh & ! function + , bl2h & ! function + , h2bl & ! function , ed_biomass & ! function , area_indices ! ! subroutine use consts_coms , only : pio4 ! ! intent(in) + use ed_misc_coms , only : igrass ! ! intent(in) implicit none !----- Arguments -----------------------------------------------------------------------! @@ -648,13 +697,24 @@ subroutine update_derived_cohort_props(cpatch,ico,green_leaf_factor,lsl) real , intent(in) :: green_leaf_factor integer , intent(in) :: lsl !----- Local variables -----------------------------------------------------------------! - real :: bl - real :: bl_max + real :: bl + real :: bl_max + integer :: ipft !---------------------------------------------------------------------------------------! - - !----- Gett DBH and height from structural biomass. ------------------------------------! - cpatch%dbh(ico) = bd2dbh(cpatch%pft(ico), cpatch%bdead(ico)) - cpatch%hite(ico) = dbh2h(cpatch%pft(ico), cpatch%dbh(ico)) + + ipft = cpatch%pft(ico) + + !----- Get DBH and height --------------------------------------------------------------! + if (is_grass(ipft).and. igrass==1) then + !--Grasses get dbh_effective and height from bleaf + cpatch%dbh(ico) = bl2dbh(cpatch%bleaf(ico), ipft) + cpatch%hite(ico) = bl2h(cpatch%bleaf(ico),ipft) + else + !--Trees get dbh from bdead + cpatch%dbh(ico) = bd2dbh(ipft, cpatch%bdead(ico)) + cpatch%hite(ico) = dbh2h(ipft, cpatch%dbh(ico)) + end if + !----- Check the phenology status and whether it needs to change. ----------------------! select case (cpatch%phenology_status(ico)) @@ -667,8 +727,14 @@ subroutine update_derived_cohort_props(cpatch,ico,green_leaf_factor,lsl) cpatch%elongf(ico) = 1.0 end select - bl_max = dbh2bl(cpatch%dbh(ico),cpatch%pft(ico)) & - * green_leaf_factor * cpatch%elongf(ico) + + if (is_grass(ipft).and. igrass==1) then + !---ALS== not sure what maximum leaf should be for grasses? Use actual leaves for now. + bl_max = cpatch%bleaf(ico) + else + bl_max = dbh2bl(cpatch%dbh(ico),ipft) * green_leaf_factor * cpatch%elongf(ico) + end if + !------------------------------------------------------------------------------------! ! If LEAF biomass is not the maximum, set it to 1 (leaves partially flushed), ! ! otherwise, set it to 0 (leaves are fully flushed). ! @@ -685,16 +751,15 @@ subroutine update_derived_cohort_props(cpatch,ico,green_leaf_factor,lsl) call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico), cpatch%hite(ico),cpatch%pft(ico) & ,cpatch%sla(ico),cpatch%lai(ico),cpatch%wai(ico),cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) !----- Finding the new basal area and above-ground biomass. ----------------------------! - cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico),cpatch%bleaf(ico) & - ,cpatch%pft(ico),cpatch%hite(ico) ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%basarea(ico)= pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !----- Update rooting depth ------------------------------------------------------------! - cpatch%krdepth(ico) = dbh2krdepth(cpatch%hite(ico),cpatch%dbh(ico),cpatch%pft(ico),lsl) + cpatch%krdepth(ico) = dbh2krdepth(cpatch%hite(ico),cpatch%dbh(ico),ipft,lsl) return end subroutine update_derived_cohort_props @@ -760,10 +825,8 @@ subroutine update_vital_rates(cpatch,ico,ilu,dbh_in,bdead_in,balive_in,hite_in,b !----- Find the new basal area and above-ground biomass. -------------------------------! cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico) ) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) !---------------------------------------------------------------------------------------! ! Change the agb growth to kgC/plant/year, basal area to cm2/plant/year, and DBH ! diff --git a/ED/src/init/ed_init_atm.F90 b/ED/src/init/ed_init_atm.F90 index 5d87f52e6..c9576a7d7 100644 --- a/ED/src/init/ed_init_atm.F90 +++ b/ED/src/init/ed_init_atm.F90 @@ -168,7 +168,7 @@ subroutine ed_init_atm() call calc_veg_hcap( cpatch%bleaf (ico) , cpatch%bdead (ico) & - , cpatch%bsapwood (ico) , cpatch%nplant (ico) & + , cpatch%bsapwooda (ico) , cpatch%nplant (ico) & , cpatch%pft (ico) , cpatch%leaf_hcap(ico) & , cpatch%wood_hcap (ico) ) diff --git a/ED/src/init/ed_nbg_init.f90 b/ED/src/init/ed_nbg_init.f90 index 281fea13e..618bfa947 100644 --- a/ED/src/init/ed_nbg_init.f90 +++ b/ED/src/init/ed_nbg_init.f90 @@ -115,14 +115,15 @@ subroutine init_nbg_cohorts(csite,lsl,ipa_a,ipa_z) , include_pft & ! intent(in) , include_these_pft & ! intent(in) , include_pft_ag & ! intent(in) - , init_density ! ! intent(in) + , init_density & ! intent(in) + , agf_bs ! ! intent(in) use consts_coms , only : t3ple & ! intent(in) , pio4 & ! intent(in) , kgom2_2_tonoha & ! intent(in) , tonoha_2_kgom2 ! ! intent(in) use allometry , only : h2dbh & ! function , dbh2bd & ! function - , dbh2bl & ! function + , size2bl & ! function , ed_biomass & ! function , area_indices ! ! subroutine use fuse_fiss_utils , only : sort_cohorts ! ! subroutine @@ -207,7 +208,7 @@ subroutine init_nbg_cohorts(csite,lsl,ipa_a,ipa_z) cpatch%phenology_status(ico) = 0 cpatch%dbh(ico) = h2dbh(cpatch%hite(ico),ipft) cpatch%bdead(ico) = dbh2bd(cpatch%dbh(ico),ipft) - cpatch%bleaf(ico) = dbh2bl(cpatch%dbh(ico),ipft) + cpatch%bleaf(ico) = size2bl(cpatch%dbh(ico),cpatch%hite(ico),ipft) cpatch%sla(ico) = sla(ipft) @@ -216,22 +217,22 @@ subroutine init_nbg_cohorts(csite,lsl,ipa_a,ipa_z) cpatch%balive(ico) = cpatch%bleaf(ico) * salloc cpatch%broot(ico) = q(ipft) * cpatch%balive(ico) * salloci - cpatch%bsapwood(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & - * salloci + cpatch%bsapwooda(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & + * salloci * agf_bs(ipft) + cpatch%bsapwoodb(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & + * salloci * (1.-agf_bs(ipft)) cpatch%bstorage(ico) = 0.5 * ( cpatch%bleaf(ico) + cpatch%broot(ico) & - + cpatch%bsapwood(ico)) + + cpatch%bsapwooda(ico)+ cpatch%bsapwoodb(ico)) !----- Find the initial area indices (LAI, WAI, CAI). ----------------------------! call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico), cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & - ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Find the above-ground biomass and basal area. -----------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico),cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) cpatch%basarea(ico) = pio4 * cpatch%dbh(ico)*cpatch%dbh(ico) !----- Initialize other cohort-level variables. ----------------------------------! @@ -277,14 +278,15 @@ subroutine init_cohorts_by_layers(csite,lsl,ipa_a,ipa_z) , include_pft & ! intent(in) , include_these_pft & ! intent(in) , include_pft_ag & ! intent(in) - , init_density ! ! intent(in) + , init_density & ! intent(in) + , agf_bs ! ! intent(in) use consts_coms , only : t3ple & ! intent(in) , pio4 & ! intent(in) , kgom2_2_tonoha & ! intent(in) , tonoha_2_kgom2 ! ! intent(in) use allometry , only : h2dbh & ! function , dbh2bd & ! function - , dbh2bl & ! function + , size2bl & ! function , ed_biomass & ! function , area_indices ! ! subroutine use fuse_fiss_utils , only : sort_cohorts ! ! subroutine @@ -345,7 +347,7 @@ subroutine init_cohorts_by_layers(csite,lsl,ipa_a,ipa_z) cpatch%bstorage(ico) = 0.0 cpatch%dbh(ico) = h2dbh(cpatch%hite(ico),ipft) cpatch%bdead(ico) = dbh2bd(cpatch%dbh(ico),ipft) - cpatch%bleaf(ico) = dbh2bl(cpatch%dbh(ico),ipft) + cpatch%bleaf(ico) = size2bl(cpatch%dbh(ico),cpatch%hite(ico),ipft) cpatch%sla(ico) = sla(ipft) @@ -354,8 +356,10 @@ subroutine init_cohorts_by_layers(csite,lsl,ipa_a,ipa_z) cpatch%balive(ico) = cpatch%bleaf(ico) * salloc cpatch%broot(ico) = q(ipft) * cpatch%balive(ico) * salloci - cpatch%bsapwood(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & - * salloci + cpatch%bsapwooda(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & + * salloci * agf_bs(ipft) + cpatch%bsapwoodb(ico) = qsw(ipft) * cpatch%hite(ico) * cpatch%balive(ico) & + * salloci * (1.-agf_bs(ipft)) !----- NPlant is defined such that the cohort LAI is equal to LAI0 cpatch%nplant(ico) = lai0 / (cpatch%bleaf(ico) * cpatch%sla(ico)) @@ -364,13 +368,11 @@ subroutine init_cohorts_by_layers(csite,lsl,ipa_a,ipa_z) call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico), cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & - ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%wai(ico),cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Find the above-ground biomass and basal area. -----------------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico),cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) !----- Initialize other cohort-level variables. ----------------------------------! @@ -415,6 +417,7 @@ subroutine near_bare_ground_big_leaf_init(cgrid) use pft_coms , only : q & ! intent(in) , qsw & ! intent(in) , sla & ! intent(in) + , agf_bs & ! intent(in) , hgt_min & ! intent(in) , include_pft & ! intent(in) , include_these_pft & ! intent(in) @@ -534,22 +537,22 @@ subroutine near_bare_ground_big_leaf_init(cgrid) salloci = 1. / salloc cpatch%balive(ico) = cpatch%bleaf(ico) * salloc - cpatch%broot(ico) = q(ipft) * cpatch%balive(ico) * salloci - cpatch%bsapwood(ico) = qsw(ipft) * cpatch%hite(ico) & - * cpatch%balive(ico) * salloci + cpatch%broot(ico) = q(ipft) * cpatch%balive(ico) *salloci + cpatch%bsapwooda(ico) = qsw(ipft) * cpatch%hite(ico) & + * cpatch%balive(ico) * salloci *agf_bs(ipft) + cpatch%bsapwoodb(ico) = qsw(ipft) * cpatch%hite(ico) & + * cpatch%balive(ico) * salloci *(1.-agf_bs(ipft)) !----- Find the initial area indices (LAI, WAI, CAI). ----------------------! call area_indices(cpatch%nplant(ico),cpatch%bleaf(ico),cpatch%bdead(ico) & ,cpatch%balive(ico),cpatch%dbh(ico), cpatch%hite(ico) & ,cpatch%pft(ico),cpatch%sla(ico),cpatch%lai(ico) & ,cpatch%wai(ico),cpatch%crown_area(ico) & - ,cpatch%bsapwood(ico)) + ,cpatch%bsapwooda(ico)) !----- Find the above-ground biomass and basal area. -----------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico),cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) cpatch%basarea(ico) = pio4 * cpatch%dbh(ico)*cpatch%dbh(ico) !----- Initialize other cohort-level variables. ----------------------------! diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index ca72a7b93..071168f51 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -1578,8 +1578,8 @@ subroutine init_pft_photo_params() dark_respiration_factor(11) = gamma_c3 dark_respiration_factor(12) = gamma_c3 dark_respiration_factor(13) = gamma_c3 - dark_respiration_factor(14) = gamma_c3 - dark_respiration_factor(15) = gamma_c3 + dark_respiration_factor(14) = gamma_c3 ! why are these not c4? + dark_respiration_factor(15) = gamma_c3 ! why are these not c4? dark_respiration_factor(16) = gamma_c3 dark_respiration_factor(17) = gamma_c3 * 1.2 !---------------------------------------------------------------------------------------! @@ -2197,6 +2197,7 @@ subroutine init_pft_alloc_params() use ed_max_dims , only : n_pft & ! intent(in) , str_len ! ! intent(in) use ed_misc_coms , only : iallom & ! intent(in) + , igrass & ! intent(in) , ibigleaf ! ! intent(in) use detailed_coms, only : idetailed ! ! intent(in) implicit none @@ -2204,6 +2205,7 @@ subroutine init_pft_alloc_params() integer :: ipft integer :: n real :: aux + real :: init_density_grass real :: init_bleaf logical :: write_allom !----- Constants shared by both bdead and bleaf (tropical PFTs) ------------------------! @@ -2320,7 +2322,7 @@ subroutine init_pft_alloc_params() sla_slope = -0.46 !----- [KIM] - new tropical parameters. ------------------------------------------------! - SLA( 1) = 21.0 ! 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate( 1))) * sla_scale + SLA( 1) = 22.7 !--value from Mike Dietze: mean: 22.7, median 19.1, 95% CI: 5.7, 78.6 SLA( 2) = 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate( 2))) * sla_scale SLA( 3) = 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate( 3))) * sla_scale SLA( 4) = 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate( 4))) * sla_scale @@ -2335,7 +2337,7 @@ subroutine init_pft_alloc_params() SLA(13) = 22.0 SLA(14) = 21.0 ! 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate(14))) * sla_scale SLA(15) = 21.0 ! 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate(15))) * sla_scale - SLA(16) = 21.0 ! 10.0**(sla_inter + sla_slope * log10(12.0/leaf_turnover_rate(16))) * sla_scale + SLA(16) = 22.7 !--value from Mike Dietze: mean: 22.7, median 19.1, 95% CI: 5.7, 78.6 SLA(17) = 10.0 !---------------------------------------------------------------------------------------! @@ -2394,6 +2396,28 @@ subroutine init_pft_alloc_params() + !---------------------------------------------------------------------------------------! + ! Initial density of plants, for near-bare-ground simulations [# of individuals/m2] ! + !---------------------------------------------------------------------------------------! + if (igrass==1) then + init_density_grass = 1. + else + init_density_grass = 0.1 + end if + + init_density(1) = init_density_grass + init_density(2:4) = 0.1 + init_density(5) = 0.1 + init_density(6:8) = 0.1 + init_density(9:11) = 0.1 + init_density(12:13) = 0.1 + init_density(14:15) = 0.1 + init_density(16) = init_density_grass + init_density(17) = 0.1 + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! ! DBH/height allometry parameters. ! ! ! @@ -2954,13 +2978,13 @@ subroutine init_pft_leaf_params() phenology(16) = 1 phenology(17) = 0 case (2) - phenology(1) = 4 + phenology(1) = 0 phenology(2:4) = 4 - phenology(5) = 4 + phenology(5) = 0 phenology(6:8) = 0 phenology(9:11) = 2 phenology(12:15) = 4 - phenology(16) = 4 + phenology(16) = 0 phenology(17) = 0 case (3) phenology(1) = 4 @@ -3125,6 +3149,7 @@ subroutine init_pft_derived_params() use allometry , only : h2dbh & ! function , dbh2h & ! function , dbh2bl & ! function + , size2bl & ! function , dbh2bd ! ! function use ed_therm_lib , only : calc_veg_hcap ! ! function implicit none @@ -3182,7 +3207,7 @@ subroutine init_pft_derived_params() !----- Find the DBH and carbon pools associated with a newly formed recruit. --------! dbh = h2dbh(hgt_min(ipft),ipft) - bleaf_min = dbh2bl(dbh,ipft) + bleaf_min = size2bl(dbh,hgt_min(ipft),ipft) broot_min = bleaf_min * q(ipft) bsapwood_min = bleaf_min * qsw(ipft) * hgt_min(ipft) balive_min = bleaf_min + broot_min + bsapwood_min @@ -3220,6 +3245,12 @@ subroutine init_pft_derived_params() min_recruit_size(ipft) = min_plant_dens * one_plant_c(ipft) !------------------------------------------------------------------------------------! + write (unit=*,fmt='(a,1x,es12.4)') ' - min_recruit_size: ',min_recruit_size(ipft) + write (unit=*,fmt='(a,1x,es12.4)') ' - min_plant_den: ',min_plant_dens + write (unit=*,fmt='(a,1x,es12.4)') ' - balive_min: ',balive_min + write (unit=*,fmt='(a,1x,es12.4)') ' - bdead_min: ',bdead_min + write (unit=*,fmt='(a,1x,es12.4)') ' - dbh: ',dbh + write (unit=*,fmt='(a,1x,es12.4)') ' - bleaf_min: ',bleaf_min !------------------------------------------------------------------------------------! diff --git a/ED/src/io/average_utils.f90 b/ED/src/io/average_utils.f90 index ff905ca9c..920cd4981 100644 --- a/ED/src/io/average_utils.f90 +++ b/ED/src/io/average_utils.f90 @@ -444,14 +444,15 @@ subroutine reset_averaged_vars(cgrid) cgrid%avg_balive (ipy) = 0.0 cgrid%avg_bleaf (ipy) = 0.0 cgrid%avg_broot (ipy) = 0.0 - cgrid%avg_bsapwood (ipy) = 0.0 - cgrid%avg_bstorage (ipy) = 0.0 + cgrid%avg_bsapwooda (ipy) = 0.0 + cgrid%avg_bsapwoodb (ipy) = 0.0 + cgrid%avg_bstorage (ipy) = 0.0 cgrid%avg_bseeds (ipy) = 0.0 cgrid%avg_fsc (ipy) = 0.0 - cgrid%avg_ssc (ipy) = 0.0 + cgrid%avg_ssc (ipy) = 0.0 cgrid%avg_stsc (ipy) = 0.0 cgrid%avg_fsn (ipy) = 0.0 - cgrid%avg_msn (ipy) = 0.0 + cgrid%avg_msn (ipy) = 0.0 cpoly => cgrid%polygon(ipy) @@ -4057,7 +4058,6 @@ subroutine update_ed_yearly_vars(cgrid) use ed_state_vars,only:edtype,polygontype,sitetype,patchtype use ed_max_dims, only: n_pft, n_dbh use consts_coms, only: pi1 - use allometry, only: ed_biomass implicit none diff --git a/ED/src/io/ed_init_full_history.F90 b/ED/src/io/ed_init_full_history.F90 index 3c29839db..c2bc9f3e2 100644 --- a/ED/src/io/ed_init_full_history.F90 +++ b/ED/src/io/ed_init_full_history.F90 @@ -2556,7 +2556,8 @@ end subroutine hdf_getslab_i call hdf_getslab_i(cpatch%phenology_status,'PHENOLOGY_STATUS ',dsetrank,iparallel,.true.) call hdf_getslab_r(cpatch%balive,'BALIVE ',dsetrank,iparallel,.true.) call hdf_getslab_r(cpatch%broot,'BROOT ',dsetrank,iparallel,.true.) - call hdf_getslab_r(cpatch%bsapwood,'BSAPWOOD ',dsetrank,iparallel,.true.) + call hdf_getslab_r(cpatch%bsapwooda,'BSAPWOODA ',dsetrank,iparallel,.true.) + call hdf_getslab_r(cpatch%bsapwoodb,'BSAPWOODB ',dsetrank,iparallel,.true.) call hdf_getslab_r(cpatch%lai,'LAI_CO ',dsetrank,iparallel,.true.) call hdf_getslab_r(cpatch%llspan,'LLSPAN ',dsetrank,iparallel,.true.) @@ -2573,7 +2574,7 @@ end subroutine hdf_getslab_i if (all(cpatch%crown_area == 0.)) then do ico= 1,cpatch%ncohorts cpatch%crown_area(ico) = min(1.0, cpatch%nplant(ico) * dbh2ca(cpatch%dbh(ico) & - ,cpatch%sla(ico),cpatch%pft(ico))) + ,cpatch%hite(ico),cpatch%sla(ico),cpatch%pft(ico))) end do end if diff --git a/ED/src/io/ed_load_namelist.f90 b/ED/src/io/ed_load_namelist.f90 index 7a577b0da..7bc8f07ab 100644 --- a/ED/src/io/ed_load_namelist.f90 +++ b/ED/src/io/ed_load_namelist.f90 @@ -116,6 +116,7 @@ subroutine copy_nl(copy_type) , lwidth_nltree & ! intent(out) , q10_c3 & ! intent(out) , q10_c4 & ! intent(out) + , lturnover_grass & ! intent(out) , quantum_efficiency_T ! ! intent(out) use phenology_coms , only : iphen_scheme & ! intent(out) , iphenys1 & ! intent(out) @@ -193,6 +194,7 @@ subroutine copy_nl(copy_type) , unitstate & ! intent(out) , event_file & ! intent(out) , iallom & ! intent(out) + , igrass & ! intent(out) , min_site_area ! ! intent(out) use grid_coms , only : time & ! intent(out) , centlon & ! intent(out) @@ -342,6 +344,7 @@ subroutine copy_nl(copy_type) ibranch_thermo = nl%ibranch_thermo iphysiol = nl%iphysiol iallom = nl%iallom + igrass = nl%igrass iphen_scheme = nl%iphen_scheme repro_scheme = nl%repro_scheme lapse_scheme = nl%lapse_scheme @@ -380,6 +383,7 @@ subroutine copy_nl(copy_type) lwidth_nltree = nl%lwidth_nltree q10_c3 = nl%q10_c3 q10_c4 = nl%q10_c4 + lturnover_grass = nl%lturnover_grass thetacrit = nl%thetacrit quantum_efficiency_T = nl%quantum_efficiency_T radint = nl%radint diff --git a/ED/src/io/ed_opspec.F90 b/ED/src/io/ed_opspec.F90 index 53110f95c..019f08da5 100644 --- a/ED/src/io/ed_opspec.F90 +++ b/ED/src/io/ed_opspec.F90 @@ -874,7 +874,7 @@ subroutine ed_opspec_times ! This is fine but now outstate must be set exactly as frqstate. If the user wasn't ! ! aware of this, print an informative banner. ! !------------------------------------------------------------------------------------! - elseif (outstate /= 0. .or. outstate > frqstate) then + elseif (outstate /= 0. .and. outstate > frqstate) then outstate = frqstate nrec_state = 1 write (unit=*,fmt='(a)') ' ' @@ -1115,6 +1115,7 @@ subroutine ed_opspec_misc , ibigleaf & ! intent(in) , integration_scheme & ! intent(in) , iallom & ! intent(in) + , igrass & ! intent(in) , min_site_area ! ! intent(in) use canopy_air_coms , only : icanturb & ! intent(in) , isfclyrm & ! intent(in) @@ -1552,6 +1553,14 @@ subroutine ed_opspec_misc ifaterr = ifaterr +1 end if + if (igrass < 0 .or. igrass > 1) then + write (reason,fmt='(a,1x,i4,a)') & + 'Invalid IGRASS, it must be between 0 and 1. Yours is set to' & + ,igrass,'...' + call opspec_fatal(reason,'opspec_misc') + ifaterr = ifaterr +1 + end if + if (iphen_scheme < -1 .or. iphen_scheme > 3) then write (reason,fmt='(a,1x,i4,a)') & 'Invalid IPHEN_SCHEME, it must be between -1 and 3. Yours is set to' & @@ -1782,7 +1791,7 @@ subroutine ed_opspec_misc call opspec_fatal(reason,'opspec_misc') ifaterr = ifaterr +1 end if - + if (thetacrit < -1.49 .or. thetacrit > 1.) then write (reason,fmt='(a,1x,es12.5,a)') & 'Invalid THETACRIT, it must be between -1.49 and 1. Yours is set to' & diff --git a/ED/src/io/ed_read_ed10_20_history.f90 b/ED/src/io/ed_read_ed10_20_history.f90 index f01bc5a65..d12ab20f0 100644 --- a/ED/src/io/ed_read_ed10_20_history.f90 +++ b/ED/src/io/ed_read_ed10_20_history.f90 @@ -24,6 +24,7 @@ subroutine read_ed10_ed20_history_file , include_pft & ! intent(in) , include_pft_ag & ! intent(in) , pft_1st_check & ! intent(in) + , agf_bs & ! intent(in) , include_these_pft ! ! intent(in) use ed_misc_coms , only : sfilin & ! intent(in) , ied_init_mode ! ! intent(in) @@ -46,6 +47,7 @@ subroutine read_ed10_ed20_history_file , h2dbh & ! function , dbh2bd & ! function , dbh2bl & ! function + , size2bl & ! function , ed_biomass & ! function , area_indices ! ! subroutine use fuse_fiss_utils, only : sort_cohorts & ! subroutine @@ -693,13 +695,17 @@ subroutine read_ed10_ed20_history_file ! Use allometry to define leaf and the other live biomass ! ! pools. ! !------------------------------------------------------------------! - cpatch%bleaf(ic2) = dbh2bl(dbh(ic),ipft(ic)) - cpatch%balive(ic2) = cpatch%bleaf(ic2) * (1.0 + q(ipft(ic)) & - + qsw(ipft(ic)) * cpatch%hite(ic2)) - cpatch%broot(ic2) = cpatch%balive(ic2) * q(ipft(ic)) & - / ( 1.0 + q(ipft(ic)) + qsw(ipft(ic)) & - * cpatch%hite(ic2)) - cpatch%bsapwood(ic2) = cpatch%balive(ic2) & + cpatch%bleaf(ic2) = size2bl(dbh(ic), hite(ic),ipft(ic)) + cpatch%balive(ic2) = cpatch%bleaf(ic2) * (1.0 + q(ipft(ic)) & + + qsw(ipft(ic)) * cpatch%hite(ic2)) + cpatch%broot(ic2) = cpatch%balive(ic2) * q(ipft(ic)) & + / ( 1.0 + q(ipft(ic)) + qsw(ipft(ic)) & + * cpatch%hite(ic2)) + cpatch%bsapwooda(ic2) = agf_bs(ipft(ic)) * cpatch%balive(ic2) & + * qsw(ipft(ic))* cpatch%hite(ic2) & + / ( 1.0 + q(ipft(ic)) + qsw(ipft(ic)) & + * cpatch%hite(ic2)) + cpatch%bsapwoodb(ic2) = (1.-agf_bs(ipft(ic))) * cpatch%balive(ic2) & * qsw(ipft(ic))* cpatch%hite(ic2) & / ( 1.0 + q(ipft(ic)) + qsw(ipft(ic)) & * cpatch%hite(ic2)) @@ -721,7 +727,7 @@ subroutine read_ed10_ed20_history_file ,cpatch%dbh(ic2), cpatch%hite(ic2) & ,cpatch%pft(ic2), SLA(cpatch%pft(ic2)) & ,cpatch%lai(ic2), cpatch%wai(ic2) & - ,cpatch%crown_area(ic2),cpatch%bsapwood(ic2)) + ,cpatch%crown_area(ic2),cpatch%bsapwooda(ic2)) !----- Initialise the carbon balance. -----------------------------! cpatch%cb (1:12,ic2) = cb(1:12,ic) @@ -730,10 +736,8 @@ subroutine read_ed10_ed20_history_file cpatch%cb_max( 13,ic2) = 0.0 !----- Above ground biomass, use the allometry. -------------------! - cpatch%agb(ic2) = ed_biomass(cpatch%bdead(ic2),cpatch%balive(ic2) & - ,cpatch%bleaf(ic2),cpatch%pft(ic2) & - ,cpatch%hite(ic2),cpatch%bstorage(ic2) & - ,cpatch%bsapwood(ic2)) + cpatch%agb(ic2) = ed_biomass(cpatch%bdead(ic2),cpatch%bleaf(ic2) & + ,cpatch%bsapwooda(ic2),cpatch%pft(ic2)) cpatch%basarea(ic2) = pio4 * cpatch%dbh(ic2) * cpatch%dbh(ic2) !----- Growth rates, start with zero. -----------------------------! diff --git a/ED/src/io/ed_read_ed21_history.F90 b/ED/src/io/ed_read_ed21_history.F90 index 7af64da5a..cf31d5671 100644 --- a/ED/src/io/ed_read_ed21_history.F90 +++ b/ED/src/io/ed_read_ed21_history.F90 @@ -22,6 +22,7 @@ subroutine read_ed21_history_file , include_pft & ! intent(in) , include_pft_ag & ! intent(in) , pft_1st_check & ! intent(in) + , agf_bs & ! intent(in) , include_these_pft ! ! intent(in) use ed_misc_coms , only : sfilin & ! intent(in) , current_time & ! intent(in) @@ -53,6 +54,7 @@ subroutine read_ed21_history_file use allometry , only : area_indices & ! function , ed_biomass & ! function , bd2dbh & ! function + , size2bl & ! function , dbh2h & ! function , dbh2bd & ! function , dbh2bl ! ! function @@ -596,15 +598,19 @@ subroutine read_ed21_history_file cpatch%bdead(ico) = dbh2bd(cpatch%dbh(ico),cpatch%pft (ico)) end if - cpatch%bleaf(ico) = dbh2bl(cpatch%dbh(ico),cpatch%pft(ico)) + cpatch%bleaf(ico) = size2bl(cpatch%dbh(ico),cpatch%hite(ico) & + ,cpatch%pft(ico)) !----- Find the other pools. -----------------------------------! salloc = (1.0 + q(ipft) + qsw(ipft) * cpatch%hite(ico)) salloci = 1.0 / salloc cpatch%balive (ico) = cpatch%bleaf(ico) * salloc cpatch%broot (ico) = cpatch%balive(ico) * q(ipft) * salloci - cpatch%bsapwood(ico) = cpatch%balive(ico) * qsw(ipft) & - * cpatch%hite(ico) * salloci + cpatch%bsapwooda(ico) = cpatch%balive(ico) * qsw(ipft) & + * cpatch%hite(ico) * salloci * agf_bs(ipft) + cpatch%bsapwoodb(ico) = cpatch%balive(ico) * qsw(ipft) & + * cpatch%hite(ico) * salloci & + * (1.-agf_bs(ipft)) cpatch%bstorage(ico) = 0.0 cpatch%phenology_status(ico) = 0 end do @@ -695,24 +701,28 @@ subroutine read_ed21_history_file cpatch%balive(ico) < tiny_biomass) then cpatch%balive(ico) = tiny_biomass end if - if (cpatch%bleaf(ico) > 0. .and. & - cpatch%bleaf(ico) < tiny_biomass) then + if (cpatch%bleaf(ico) > 0. .and. & + cpatch%bleaf(ico) < tiny_biomass) then cpatch%bleaf(ico) = tiny_biomass end if - if (cpatch%broot(ico) > 0. .and. & - cpatch%broot(ico) < tiny_biomass) then + if (cpatch%broot(ico) > 0. .and. & + cpatch%broot(ico) < tiny_biomass) then cpatch%broot(ico) = tiny_biomass end if - if (cpatch%bsapwood(ico) > 0. .and. & - cpatch%bsapwood(ico) < tiny_biomass) then - cpatch%bsapwood(ico) = tiny_biomass + if (cpatch%bsapwooda(ico) > 0. .and. & + cpatch%bsapwooda(ico) < tiny_biomass) then + cpatch%bsapwooda(ico) = tiny_biomass end if - if (cpatch%bdead(ico) > 0. .and. & - cpatch%bdead(ico) < tiny_biomass) then + if (cpatch%bsapwoodb(ico) > 0. .and. & + cpatch%bsapwoodb(ico) < tiny_biomass) then + cpatch%bsapwoodb(ico) = tiny_biomass + end if + if (cpatch%bdead(ico) > 0. .and. & + cpatch%bdead(ico) < tiny_biomass) then cpatch%bdead(ico) = tiny_biomass end if - if (cpatch%bstorage(ico) > 0. .and. & - cpatch%bstorage(ico) < tiny_biomass) then + if (cpatch%bstorage(ico) > 0. .and. & + cpatch%bstorage(ico) < tiny_biomass) then cpatch%bstorage(ico) = tiny_biomass end if !---------------------------------------------------------------! @@ -721,14 +731,8 @@ subroutine read_ed21_history_file !----- Compute the above-ground biomass. -----------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico) & - ,cpatch%balive(ico) & - ,cpatch%bleaf(ico) & - ,cpatch%pft(ico) & - ,cpatch%hite(ico) & - ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico) ) - + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico)& + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) @@ -738,7 +742,7 @@ subroutine read_ed21_history_file ,cpatch%dbh(ico), cpatch%hite(ico) & ,cpatch%pft(ico), SLA(cpatch%pft(ico)) & ,cpatch%lai(ico),cpatch%wai(ico) & - ,cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Update the derived patch-level variables. ---------------! @@ -859,6 +863,7 @@ subroutine read_ed21_history_unstruct , include_pft_ag & ! intent(in) , pft_1st_check & ! intent(in) , include_these_pft & ! intent(in) + , agf_bs & ! intent(in) , min_cohort_size ! ! intent(in) use ed_misc_coms , only : sfilin & ! intent(in) , current_time & ! intent(in) @@ -1752,11 +1757,14 @@ subroutine read_ed21_history_unstruct !----- Find the other pools. -----------------------------------! salloc = (1.0 + q(ipft) + qsw(ipft) * cpatch%hite(ico)) salloci = 1.0 / salloc - cpatch%balive (ico) = cpatch%bleaf(ico) * salloc - cpatch%broot (ico) = cpatch%balive(ico) * q(ipft) * salloci - cpatch%bsapwood(ico) = cpatch%balive(ico) * qsw(ipft) & - * cpatch%hite(ico) * salloci - cpatch%bstorage(ico) = 0.0 + cpatch%balive (ico) = cpatch%bleaf(ico) * salloc + cpatch%broot (ico) = cpatch%balive(ico) * q(ipft) * salloci + cpatch%bsapwooda(ico) = cpatch%balive(ico) * qsw(ipft) & + * cpatch%hite(ico) * salloci * agf_bs(ipft) + cpatch%bsapwoodb(ico) = cpatch%balive(ico) * qsw(ipft) & + * cpatch%hite(ico) * salloci & + * (1.-agf_bs(ipft)) + cpatch%bstorage(ico) = 0.0 cpatch%phenology_status(ico) = 0 end do @@ -1854,9 +1862,13 @@ subroutine read_ed21_history_unstruct cpatch%broot(ico) < tiny_biomass) then cpatch%broot(ico) = tiny_biomass end if - if (cpatch%bsapwood(ico) > 0. .and. & - cpatch%bsapwood(ico) < tiny_biomass) then - cpatch%bsapwood(ico) = tiny_biomass + if (cpatch%bsapwooda(ico) > 0. .and. & + cpatch%bsapwooda(ico) < tiny_biomass) then + cpatch%bsapwooda(ico) = tiny_biomass + end if + if (cpatch%bsapwoodb(ico) > 0. .and. & + cpatch%bsapwoodb(ico) < tiny_biomass) then + cpatch%bsapwoodb(ico) = tiny_biomass end if if (cpatch%bdead(ico) > 0. .and. & cpatch%bdead(ico) < tiny_biomass) then @@ -1870,12 +1882,8 @@ subroutine read_ed21_history_unstruct !----- Compute the above-ground biomass. -----------------------! - cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico) & - ,cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico) & - ,cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) + cpatch%agb(ico) = ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico)& + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) cpatch%basarea(ico) = pio4 * cpatch%dbh(ico) * cpatch%dbh(ico) @@ -1885,7 +1893,7 @@ subroutine read_ed21_history_unstruct ,cpatch%dbh(ico),cpatch%hite(ico) & ,cpatch%pft(ico),SLA(cpatch%pft(ico)) & ,cpatch%lai(ico),cpatch%wai(ico) & - ,cpatch%crown_area(ico),cpatch%bsapwood(ico)) + ,cpatch%crown_area(ico),cpatch%bsapwooda(ico)) !----- Update the derived patch-level variables. ---------------! diff --git a/ED/src/io/edio.f90 b/ED/src/io/edio.f90 index 4308119a2..243e455ca 100644 --- a/ED/src/io/edio.f90 +++ b/ED/src/io/edio.f90 @@ -393,7 +393,8 @@ subroutine spatial_averages cgrid%avg_balive (ipy) = 0.0 cgrid%avg_bleaf (ipy) = 0.0 cgrid%avg_broot (ipy) = 0.0 - cgrid%avg_bsapwood (ipy) = 0.0 + cgrid%avg_bsapwooda (ipy) = 0.0 + cgrid%avg_bsapwoodb (ipy) = 0.0 cgrid%avg_bdead (ipy) = 0.0 cgrid%avg_bstorage (ipy) = 0.0 cgrid%avg_bseeds (ipy) = 0.0 @@ -810,8 +811,13 @@ subroutine spatial_averages * csite%area(ipa)*cpoly%area(isi) & * site_area_i * poly_area_i - cgrid%avg_bsapwood(ipy) = cgrid%avg_bsapwood(ipy) & - + sum(cpatch%bsapwood*cpatch%nplant) & + cgrid%avg_bsapwooda(ipy) = cgrid%avg_bsapwooda(ipy) & + + sum(cpatch%bsapwooda*cpatch%nplant) & + * csite%area(ipa)*cpoly%area(isi) & + * site_area_i * poly_area_i + + cgrid%avg_bsapwoodb(ipy) = cgrid%avg_bsapwoodb(ipy) & + + sum(cpatch%bsapwoodb*cpatch%nplant) & * csite%area(ipa)*cpoly%area(isi) & * site_area_i * poly_area_i diff --git a/ED/src/memory/ed_misc_coms.f90 b/ED/src/memory/ed_misc_coms.f90 index 81df19b57..a4893d611 100644 --- a/ED/src/memory/ed_misc_coms.f90 +++ b/ED/src/memory/ed_misc_coms.f90 @@ -218,4 +218,11 @@ Module ed_misc_coms ! the height is 35.0m. !---------------------------------------------------------------------------------------! + + !----- Namelist option for the new grass scheme. (ALS, Feb 2012) -----------------------! + integer :: igrass ! 0 -- Original ED-2.1 grasses (aka bonzai grass) + ! 1 -- New grass scheme which has bdead = 0, height = f(bleaf), and + ! growth occurs daily + !---------------------------------------------------------------------------------------! + end module ed_misc_coms diff --git a/ED/src/memory/ed_state_vars.f90 b/ED/src/memory/ed_state_vars.f90 index 49ed3cd1e..e314d272e 100644 --- a/ED/src/memory/ed_state_vars.f90 +++ b/ED/src/memory/ed_state_vars.f90 @@ -106,9 +106,12 @@ module ed_state_vars ! Biomass of fine roots (kgC/plant) real, pointer, dimension(:) :: broot - - ! Biomass of sapwood (kgC/plant) - real, pointer, dimension(:) :: bsapwood + + ! Biomass of sapwood above ground (kgC/plant) + real, pointer, dimension(:) :: bsapwooda + + ! Biomass of sapwood below ground (kgC/plant) + real, pointer, dimension(:) :: bsapwoodb ! Leaf area index (m2 leaf / m2 ground) real ,pointer,dimension(:) :: lai @@ -1374,7 +1377,8 @@ module ed_state_vars real,pointer,dimension(:) :: avg_balive real,pointer,dimension(:) :: avg_bleaf real,pointer,dimension(:) :: avg_broot - real,pointer,dimension(:) :: avg_bsapwood + real,pointer,dimension(:) :: avg_bsapwooda + real,pointer,dimension(:) :: avg_bsapwoodb real,pointer,dimension(:) :: avg_bdead real,pointer,dimension(:) :: avg_fsn real,pointer,dimension(:) :: avg_msn @@ -1655,7 +1659,8 @@ module ed_state_vars real,pointer,dimension(:) :: avg_balive real,pointer,dimension(:) :: avg_bleaf real,pointer,dimension(:) :: avg_broot - real,pointer,dimension(:) :: avg_bsapwood + real,pointer,dimension(:) :: avg_bsapwooda + real,pointer,dimension(:) :: avg_bsapwoodb real,pointer,dimension(:) :: avg_bseeds real,pointer,dimension(:) :: avg_bstorage real,pointer,dimension(:) :: avg_fsc @@ -2381,7 +2386,8 @@ subroutine allocate_edtype(cgrid,npolygons) allocate(cgrid%avg_balive (npolygons)) allocate(cgrid%avg_bleaf (npolygons)) allocate(cgrid%avg_broot (npolygons)) - allocate(cgrid%avg_bsapwood (npolygons)) + allocate(cgrid%avg_bsapwooda (npolygons)) + allocate(cgrid%avg_bsapwoodb (npolygons)) allocate(cgrid%avg_bstorage (npolygons)) allocate(cgrid%avg_bseeds (npolygons)) allocate(cgrid%avg_fsc (npolygons)) @@ -2859,40 +2865,40 @@ subroutine allocate_polygontype(cpoly,nsites) allocate(cpoly%avg_transp (nsites)) allocate(cpoly%avg_evap (nsites)) - allocate(cpoly%avg_smoist_gg (nzg,nsites)) - allocate(cpoly%avg_transloss (nzg,nsites)) - allocate(cpoly%avg_runoff (nsites)) - allocate(cpoly%avg_drainage (nsites)) - allocate(cpoly%avg_drainage_heat (nsites)) - allocate(cpoly%avg_rshort_gnd (nsites)) - allocate(cpoly%avg_rlong_gnd (nsites)) - allocate(cpoly%avg_ustar (nsites)) - allocate(cpoly%avg_tstar (nsites)) - allocate(cpoly%avg_qstar (nsites)) - allocate(cpoly%avg_cstar (nsites)) - allocate(cpoly%avg_carbon_ac (nsites)) - allocate(cpoly%avg_carbon_st (nsites)) - allocate(cpoly%avg_sensible_lc (nsites)) - allocate(cpoly%avg_sensible_wc (nsites)) - allocate(cpoly%avg_qwshed_vg (nsites)) - allocate(cpoly%avg_qintercepted (nsites)) - allocate(cpoly%avg_qthroughfall (nsites)) - allocate(cpoly%avg_sensible_gc (nsites)) - allocate(cpoly%avg_sensible_ac (nsites)) - allocate(cpoly%avg_sensible_gg (nzg,nsites)) + allocate(cpoly%avg_smoist_gg (nzg,nsites)) + allocate(cpoly%avg_transloss (nzg,nsites)) + allocate(cpoly%avg_runoff (nsites)) + allocate(cpoly%avg_drainage (nsites)) + allocate(cpoly%avg_drainage_heat (nsites)) + allocate(cpoly%avg_rshort_gnd (nsites)) + allocate(cpoly%avg_rlong_gnd (nsites)) + allocate(cpoly%avg_ustar (nsites)) + allocate(cpoly%avg_tstar (nsites)) + allocate(cpoly%avg_qstar (nsites)) + allocate(cpoly%avg_cstar (nsites)) + allocate(cpoly%avg_carbon_ac (nsites)) + allocate(cpoly%avg_carbon_st (nsites)) + allocate(cpoly%avg_sensible_lc (nsites)) + allocate(cpoly%avg_sensible_wc (nsites)) + allocate(cpoly%avg_qwshed_vg (nsites)) + allocate(cpoly%avg_qintercepted (nsites)) + allocate(cpoly%avg_qthroughfall (nsites)) + allocate(cpoly%avg_sensible_gc (nsites)) + allocate(cpoly%avg_sensible_ac (nsites)) + allocate(cpoly%avg_sensible_gg (nzg,nsites)) allocate(cpoly%avg_runoff_heat (nsites)) ! Fast time state diagnostics - allocate(cpoly%avg_leaf_energy (nsites)) - allocate(cpoly%avg_leaf_hcap (nsites)) - allocate(cpoly%avg_leaf_temp (nsites)) - allocate(cpoly%avg_leaf_fliq (nsites)) - allocate(cpoly%avg_leaf_water (nsites)) - allocate(cpoly%avg_wood_energy (nsites)) - allocate(cpoly%avg_wood_hcap (nsites)) - allocate(cpoly%avg_wood_temp (nsites)) - allocate(cpoly%avg_wood_fliq (nsites)) - allocate(cpoly%avg_wood_water (nsites)) + allocate(cpoly%avg_leaf_energy (nsites)) + allocate(cpoly%avg_leaf_hcap (nsites)) + allocate(cpoly%avg_leaf_temp (nsites)) + allocate(cpoly%avg_leaf_fliq (nsites)) + allocate(cpoly%avg_leaf_water (nsites)) + allocate(cpoly%avg_wood_energy (nsites)) + allocate(cpoly%avg_wood_hcap (nsites)) + allocate(cpoly%avg_wood_temp (nsites)) + allocate(cpoly%avg_wood_fliq (nsites)) + allocate(cpoly%avg_wood_water (nsites)) allocate(cpoly%avg_can_temp (nsites)) allocate(cpoly%avg_can_shv (nsites)) allocate(cpoly%avg_can_co2 (nsites)) @@ -2919,18 +2925,19 @@ subroutine allocate_polygontype(cpoly,nsites) allocate(cpoly%avg_atm_prss (nsites)) !!!NACP - allocate(cpoly%avg_sfcw_depth (nsites)) - allocate(cpoly%avg_sfcw_energy (nsites)) - allocate(cpoly%avg_sfcw_mass (nsites)) - allocate(cpoly%avg_sfcw_fracliq (nsites)) - allocate(cpoly%avg_sfcw_tempk (nsites)) + allocate(cpoly%avg_sfcw_depth (nsites)) + allocate(cpoly%avg_sfcw_energy (nsites)) + allocate(cpoly%avg_sfcw_mass (nsites)) + allocate(cpoly%avg_sfcw_fracliq (nsites)) + allocate(cpoly%avg_sfcw_tempk (nsites)) allocate(cpoly%avg_fsc (nsites)) allocate(cpoly%avg_stsc (nsites)) allocate(cpoly%avg_ssc (nsites)) allocate(cpoly%avg_balive (nsites)) allocate(cpoly%avg_bleaf (nsites)) allocate(cpoly%avg_broot (nsites)) - allocate(cpoly%avg_bsapwood (nsites)) + allocate(cpoly%avg_bsapwooda (nsites)) + allocate(cpoly%avg_bsapwoodb (nsites)) allocate(cpoly%avg_bdead (nsites)) allocate(cpoly%avg_fsn (nsites)) allocate(cpoly%avg_msn (nsites)) @@ -3263,7 +3270,8 @@ subroutine allocate_patchtype(cpatch,ncohorts) allocate(cpatch%phenology_status(ncohorts)) allocate(cpatch%balive(ncohorts)) allocate(cpatch%broot(ncohorts)) - allocate(cpatch%bsapwood(ncohorts)) + allocate(cpatch%bsapwooda(ncohorts)) + allocate(cpatch%bsapwoodb(ncohorts)) allocate(cpatch%lai(ncohorts)) allocate(cpatch%wai(ncohorts)) allocate(cpatch%crown_area(ncohorts)) @@ -3609,19 +3617,20 @@ subroutine nullify_edtype(cgrid) nullify(cgrid%avg_sfcw_mass ) nullify(cgrid%avg_sfcw_tempk ) nullify(cgrid%avg_sfcw_fracliq ) - nullify(cgrid%avg_bdead ) - nullify(cgrid%avg_balive ) + nullify(cgrid%avg_bdead ) + nullify(cgrid%avg_balive ) nullify(cgrid%avg_bleaf ) nullify(cgrid%avg_broot ) - nullify(cgrid%avg_bsapwood ) + nullify(cgrid%avg_bsapwooda ) + nullify(cgrid%avg_bsapwoodb ) nullify(cgrid%avg_bstorage ) nullify(cgrid%avg_bseeds ) - nullify(cgrid%avg_fsc ) - nullify(cgrid%avg_ssc ) - nullify(cgrid%avg_stsc ) - nullify(cgrid%avg_fsn ) - nullify(cgrid%avg_msn ) + nullify(cgrid%avg_fsc ) + nullify(cgrid%avg_ssc ) + nullify(cgrid%avg_stsc ) + nullify(cgrid%avg_fsn ) + nullify(cgrid%avg_msn ) @@ -4125,7 +4134,8 @@ subroutine nullify_polygontype(cpoly) nullify(cpoly%avg_balive ) nullify(cpoly%avg_bleaf ) nullify(cpoly%avg_broot ) - nullify(cpoly%avg_bsapwood ) + nullify(cpoly%avg_bsapwooda ) + nullify(cpoly%avg_bsapwoodb ) nullify(cpoly%avg_bdead ) nullify(cpoly%avg_fsn ) nullify(cpoly%avg_msn ) @@ -4431,7 +4441,8 @@ subroutine nullify_patchtype(cpatch) nullify(cpatch%phenology_status) nullify(cpatch%balive) nullify(cpatch%broot) - nullify(cpatch%bsapwood) + nullify(cpatch%bsapwooda) + nullify(cpatch%bsapwoodb) nullify(cpatch%lai) nullify(cpatch%wai) nullify(cpatch%crown_area) @@ -4777,7 +4788,8 @@ subroutine deallocate_edtype(cgrid) if(associated(cgrid%avg_balive )) deallocate(cgrid%avg_balive ) if(associated(cgrid%avg_broot )) deallocate(cgrid%avg_broot ) if(associated(cgrid%avg_bleaf )) deallocate(cgrid%avg_bleaf ) - if(associated(cgrid%avg_bsapwood )) deallocate(cgrid%avg_bsapwood ) + if(associated(cgrid%avg_bsapwooda )) deallocate(cgrid%avg_bsapwooda ) + if(associated(cgrid%avg_bsapwoodb )) deallocate(cgrid%avg_bsapwoodb ) if(associated(cgrid%avg_bstorage )) deallocate(cgrid%avg_bstorage ) if(associated(cgrid%avg_bseeds )) deallocate(cgrid%avg_bseeds ) if(associated(cgrid%avg_fsc )) deallocate(cgrid%avg_fsc ) @@ -5307,7 +5319,8 @@ subroutine deallocate_polygontype(cpoly) if(associated(cpoly%avg_bleaf )) deallocate(cpoly%avg_bleaf ) if(associated(cpoly%avg_broot )) deallocate(cpoly%avg_broot ) - if(associated(cpoly%avg_bsapwood )) deallocate(cpoly%avg_bsapwood ) + if(associated(cpoly%avg_bsapwooda )) deallocate(cpoly%avg_bsapwooda ) + if(associated(cpoly%avg_bsapwoodb )) deallocate(cpoly%avg_bsapwoodb ) if(associated(cpoly%avg_bstorage )) deallocate(cpoly%avg_bstorage ) if(associated(cpoly%avg_bseeds )) deallocate(cpoly%avg_bseeds ) @@ -5619,7 +5632,8 @@ subroutine deallocate_patchtype(cpatch) if(associated(cpatch%phenology_status)) deallocate(cpatch%phenology_status) if(associated(cpatch%balive)) deallocate(cpatch%balive) if(associated(cpatch%broot)) deallocate(cpatch%broot) - if(associated(cpatch%bsapwood)) deallocate(cpatch%bsapwood) + if(associated(cpatch%bsapwooda)) deallocate(cpatch%bsapwooda) + if(associated(cpatch%bsapwoodb)) deallocate(cpatch%bsapwoodb) if(associated(cpatch%lai)) deallocate(cpatch%lai) if(associated(cpatch%wai)) deallocate(cpatch%wai) if(associated(cpatch%crown_area)) deallocate(cpatch%crown_area) @@ -6468,7 +6482,8 @@ subroutine copy_patchtype_mask(patchin,patchout,mask,masksz,newsz) patchout%phenology_status(1:inc) = pack(patchin%phenology_status,mask) patchout%balive(1:inc) = pack(patchin%balive,mask) patchout%broot(1:inc) = pack(patchin%broot,mask) - patchout%bsapwood(1:inc) = pack(patchin%bsapwood,mask) + patchout%bsapwooda(1:inc) = pack(patchin%bsapwooda,mask) + patchout%bsapwoodb(1:inc) = pack(patchin%bsapwoodb,mask) patchout%lai(1:inc) = pack(patchin%lai,mask) patchout%wai(1:inc) = pack(patchin%wai,mask) patchout%crown_area(1:inc) = pack(patchin%crown_area,mask) @@ -6708,7 +6723,8 @@ subroutine copy_patchtype(patchin,patchout,ipin1,ipin2,ipout1,ipout2) patchout%phenology_status(iout) = patchin%phenology_status(iin) patchout%balive(iout) = patchin%balive(iin) patchout%broot(iout) = patchin%broot(iin) - patchout%bsapwood(iout) = patchin%bsapwood(iin) + patchout%bsapwooda(iout) = patchin%bsapwooda(iin) + patchout%bsapwoodb(iout) = patchin%bsapwoodb(iin) patchout%lai(iout) = patchin%lai(iin) patchout%wai(iout) = patchin%wai(iin) patchout%crown_area(iout) = patchin%crown_area(iin) @@ -8235,11 +8251,18 @@ subroutine filltab_edtype_p11(cgrid,igr,init,var_len,var_len_global,max_ptrs,nva call metadata_edio(nvar,igr,'Poly Avg. Biomass -- fine roots','[kgC/m2]','ipoly') end if - if (associated(cgrid%avg_bsapwood)) then + if (associated(cgrid%avg_bsapwooda)) then + nvar=nvar+1 + call vtable_edio_r(npts,cgrid%avg_bsapwooda,nvar,igr,init,cgrid%pyglob_id, & + var_len,var_len_global,max_ptrs,'AVG_BSAPWOODA :11:hist:anal:year') + call metadata_edio(nvar,igr,'Poly Avg. Biomass -- sapwood above ground','[kgC/m2]','ipoly') + end if + + if (associated(cgrid%avg_bsapwoodb)) then nvar=nvar+1 - call vtable_edio_r(npts,cgrid%avg_bsapwood,nvar,igr,init,cgrid%pyglob_id, & - var_len,var_len_global,max_ptrs,'AVG_BSAPWOOD :11:hist:anal:year') - call metadata_edio(nvar,igr,'Poly Avg. Biomass -- sapwood','[kgC/m2]','ipoly') + call vtable_edio_r(npts,cgrid%avg_bsapwoodb,nvar,igr,init,cgrid%pyglob_id, & + var_len,var_len_global,max_ptrs,'AVG_BSAPWOODB :11:hist:anal:year') + call metadata_edio(nvar,igr,'Poly Avg. Biomass -- sapwood below ground','[kgC/m2]','ipoly') end if if (associated(cgrid%avg_fsc)) then @@ -13433,10 +13456,17 @@ subroutine filltab_patchtype(igr,ipy,isi,ipa,init) call metadata_edio(nvar,igr,'No metadata available','[NA]','NA') end if - if (associated(cpatch%bsapwood)) then + if (associated(cpatch%bsapwooda)) then + nvar=nvar+1 + call vtable_edio_r(npts,cpatch%bsapwooda,nvar,igr,init,cpatch%coglob_id, & + var_len,var_len_global,max_ptrs,'BSAPWOODA :41:hist:year:dail:anal:mont:dcyc') + call metadata_edio(nvar,igr,'No metadata available','[NA]','NA') + end if + + if (associated(cpatch%bsapwoodb)) then nvar=nvar+1 - call vtable_edio_r(npts,cpatch%bsapwood,nvar,igr,init,cpatch%coglob_id, & - var_len,var_len_global,max_ptrs,'BSAPWOOD :41:hist:year:dail:mont:dcyc') + call vtable_edio_r(npts,cpatch%bsapwoodb,nvar,igr,init,cpatch%coglob_id, & + var_len,var_len_global,max_ptrs,'BSAPWOODB :41:hist:year:dail:anal:mont:dcyc') call metadata_edio(nvar,igr,'No metadata available','[NA]','NA') end if diff --git a/ED/src/memory/ename_coms.f90 b/ED/src/memory/ename_coms.f90 index 0deb6d90b..7dfa94691 100644 --- a/ED/src/memory/ename_coms.f90 +++ b/ED/src/memory/ename_coms.f90 @@ -137,6 +137,7 @@ module ename_coms integer :: ibranch_thermo integer :: iphysiol integer :: iallom + integer :: igrass integer :: iphen_scheme integer :: repro_scheme integer :: lapse_scheme @@ -175,6 +176,7 @@ module ename_coms real :: lwidth_nltree real :: q10_c3 real :: q10_c4 + real :: lturnover_grass real :: thetacrit integer :: quantum_efficiency_T integer :: n_plant_lim @@ -391,6 +393,7 @@ subroutine init_ename_vars(enl) enl%ibranch_thermo = undef_integer enl%iphysiol = undef_integer enl%iallom = undef_integer + enl%igrass = undef_integer enl%iphen_scheme = undef_integer enl%repro_scheme = undef_integer enl%lapse_scheme = undef_integer @@ -429,6 +432,7 @@ subroutine init_ename_vars(enl) enl%lwidth_nltree = undef_real enl%q10_c3 = undef_real enl%q10_c4 = undef_real + enl%lturnover_grass = undef_real enl%thetacrit = undef_real enl%quantum_efficiency_T = undef_integer enl%n_plant_lim = undef_integer diff --git a/ED/src/memory/pft_coms.f90 b/ED/src/memory/pft_coms.f90 index dc498d3bf..0b644a906 100644 --- a/ED/src/memory/pft_coms.f90 +++ b/ED/src/memory/pft_coms.f90 @@ -592,7 +592,8 @@ module pft_coms real :: bdead real :: bleaf real :: broot - real :: bsapwood + real :: bsapwooda + real :: bsapwoodb real :: balive real :: paw_avg real :: elongf @@ -633,7 +634,8 @@ subroutine zero_recruit(maxp,recruit) recruit(p)%bdead = 0. recruit(p)%bleaf = 0. recruit(p)%broot = 0. - recruit(p)%bsapwood = 0. + recruit(p)%bsapwooda = 0. + recruit(p)%bsapwoodb = 0. recruit(p)%balive = 0. recruit(p)%paw_avg = 0. recruit(p)%elongf = 0. @@ -674,7 +676,8 @@ subroutine copy_recruit(recsource,rectarget) rectarget%bdead = recsource%bdead rectarget%bleaf = recsource%bleaf rectarget%broot = recsource%broot - rectarget%bsapwood = recsource%bsapwood + rectarget%bsapwooda = recsource%bsapwooda + rectarget%bsapwoodb = recsource%bsapwoodb rectarget%balive = recsource%balive rectarget%paw_avg = recsource%paw_avg rectarget%elongf = recsource%elongf diff --git a/ED/src/memory/physiology_coms.f90 b/ED/src/memory/physiology_coms.f90 index 1edefa929..401c91c84 100644 --- a/ED/src/memory/physiology_coms.f90 +++ b/ED/src/memory/physiology_coms.f90 @@ -105,6 +105,7 @@ module physiology_coms ! only. ! ! Q10_C3 -- Q10 factor for C3 plants (used only if IPHYSIOL is set to 2 or 3). ! ! Q10_C4 -- Q10 factor for C4 plants (used only if IPHYSIOL is set to 2 or 3). ! + ! LTURNOVER_GRASS- Leave turnover rate for grasses (1/yr) ! !---------------------------------------------------------------------------------------! real(kind=4) :: vmfact_c3 real(kind=4) :: vmfact_c4 @@ -130,6 +131,7 @@ module physiology_coms real(kind=4) :: lwidth_nltree real(kind=4) :: q10_c3 real(kind=4) :: q10_c4 + real(kind=4) :: lturnover_grass !---------------------------------------------------------------------------------------! !---------------------------------------------------------------------------------------! diff --git a/ED/src/mpi/ed_mpass_init.f90 b/ED/src/mpi/ed_mpass_init.f90 index 281c48bf0..2ffe1db2b 100644 --- a/ED/src/mpi/ed_mpass_init.f90 +++ b/ED/src/mpi/ed_mpass_init.f90 @@ -132,6 +132,7 @@ subroutine ed_masterput_nl(par_run) , idateh & ! intent(in) , ndcycle & ! intent(in) , iallom & ! intent(in) + , igrass & ! intent(in) , min_site_area & ! intent(in) , attach_metadata ! ! intent(in) use canopy_air_coms , only : icanturb & ! intent(in) @@ -231,6 +232,7 @@ subroutine ed_masterput_nl(par_run) , lwidth_nltree & ! intent(in) , q10_c3 & ! intent(in) , q10_c4 & ! intent(in) + , lturnover_grass & ! intent(in) , quantum_efficiency_T ! ! intent(in) use phenology_coms , only : iphen_scheme & ! intent(in) , iphenys1 & ! intent(in) @@ -401,6 +403,7 @@ subroutine ed_masterput_nl(par_run) call MPI_Bcast(ibranch_thermo,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(iphysiol,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(iallom,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(igrass,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(iphen_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(repro_scheme,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(radint,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) @@ -442,6 +445,7 @@ subroutine ed_masterput_nl(par_run) call MPI_Bcast(lwidth_nltree,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(q10_c3,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(q10_c4,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) + call MPI_Bcast(lturnover_grass,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(thetacrit,1,MPI_REAL,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(quantum_efficiency_T,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) call MPI_Bcast(n_plant_lim,1,MPI_INTEGER,mainnum,MPI_COMM_WORLD,ierr) @@ -1221,6 +1225,7 @@ subroutine ed_nodeget_nl , idateh & ! intent(out) , ndcycle & ! intent(out) , iallom & ! intent(out) + , igrass & ! intent(out) , min_site_area & ! intent(out) , attach_metadata ! ! intent(out) use canopy_air_coms , only : icanturb & ! intent(out) @@ -1320,6 +1325,7 @@ subroutine ed_nodeget_nl , lwidth_nltree & ! intent(out) , q10_c3 & ! intent(out) , q10_c4 & ! intent(out) + , lturnover_grass & ! intent(out) , quantum_efficiency_T ! ! intent(out) use phenology_coms , only : iphen_scheme & ! intent(out) , iphenys1 & ! intent(out) @@ -1496,6 +1502,7 @@ subroutine ed_nodeget_nl call MPI_Bcast(ibranch_thermo,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(iphysiol,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(iallom,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(igrass,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(iphen_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(repro_scheme,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(radint,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) @@ -1537,6 +1544,7 @@ subroutine ed_nodeget_nl call MPI_Bcast(lwidth_nltree,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(q10_c3,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(q10_c4,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) + call MPI_Bcast(lturnover_grass,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(thetacrit,1,MPI_REAL,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(quantum_efficiency_T,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) call MPI_Bcast(n_plant_lim,1,MPI_INTEGER,master_num,MPI_COMM_WORLD,ierr) diff --git a/ED/src/utils/allometry.f90 b/ED/src/utils/allometry.f90 index 83ddda391..fcb5a0170 100644 --- a/ED/src/utils/allometry.f90 +++ b/ED/src/utils/allometry.f90 @@ -94,14 +94,18 @@ real function dbh2bd(dbh,ipft) , b2Bs_small & ! intent(in), lookup table , b1Bs_large & ! intent(in), lookup table , b2Bs_large & ! intent(in), lookup table + , is_grass & ! intent(in) , dbh_crit ! ! intent(in), lookup table + use ed_misc_coms, only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! real , intent(in) :: dbh integer, intent(in) :: ipft !------------------------------------------------------------------------------------! - if (dbh <= dbh_crit(ipft)) then + if (is_grass(ipft).and. igrass==1) then + dbh2bd = 0.0 + else if (dbh <= dbh_crit(ipft)) then dbh2bd = b1Bs_small(ipft) / C2B * dbh ** b2Bs_small(ipft) else dbh2bd = b1Bs_large(ipft) / C2B * dbh ** b2Bs_large(ipft) @@ -158,16 +162,57 @@ end function bd2dbh + !=======================================================================================! + ! Intended to replace dbh2bl and h2bl with a single generic function. This simplifies ! + ! the code in many other places for grasses. + !=======================================================================================! + real function size2bl(dbh,hite,ipft) + use pft_coms , only : dbh_crit & ! intent(in), lookup table + , C2B & ! intent(in), lookup table + , b1Bl & ! intent(in), lookup table + , b2Bl & ! intent(in), lookup table + , hgt_max & ! intent(in), lookup table + , is_grass ! ! intent(in) + use ed_misc_coms, only : igrass ! ! intent(in) + + + implicit none + !----- Arguments --------------------------------------------------------------------! + real , intent(in) :: dbh + real , intent(in) :: hite + integer, intent(in) :: ipft + !----- Local variables --------------------------------------------------------------! + real :: mdbh + !------------------------------------------------------------------------------------! + + if (is_grass(ipft) .and. igrass==1) then + !-- use height for grasses + mdbh = min(h2dbh(hite,ipft),dbh_crit(ipft)) + else + !--use dbh for trees + mdbh = min(dbh,dbh_crit(ipft)) + end if + + size2bl = b1Bl(ipft) / C2B * mdbh ** b2Bl(ipft) + + return + end function size2bl + !=======================================================================================! + !=======================================================================================! + !=======================================================================================! !=======================================================================================! ! This function determines the maximum leaf biomass (kgC/plant) real function dbh2bl(dbh,ipft) - use pft_coms , only : dbh_crit & ! intent(in), lookup table - , C2B & ! intent(in) - , b1Bl & ! intent(in), lookup table - , b2Bl ! ! intent(in), lookup table + use pft_coms , only : dbh_crit & ! intent(in), lookup table + , C2B & ! intent(in), lookup table + , b1Bl & ! intent(in), lookup table + , b2Bl & ! intent(in), lookup table + , hgt_max & ! intent(in), lookup table + , is_grass ! ! intent(in) + use ed_misc_coms, only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! @@ -176,28 +221,109 @@ real function dbh2bl(dbh,ipft) !----- Local variables --------------------------------------------------------------! real :: mdbh !------------------------------------------------------------------------------------! - - + + if (is_grass(ipft) .and. igrass==1) then + !--- The grasses really shouldent use this function at all - use size2bl + ! test grasses against maximum height rather than maximum dbh + mdbh = min(dbh, h2dbh(hgt_max(ipft),ipft)) + else + mdbh = min(dbh,dbh_crit(ipft)) + end if + !------------------------------------------------------------------------------------! - ! Make sure bleaf won't keep growing once the plant hits the maximum height. ! + ! Find maximum leaf biomass. ! !------------------------------------------------------------------------------------! - mdbh = min(dbh,dbh_crit(ipft)) + dbh2bl = b1Bl(ipft) / C2B * mdbh ** b2Bl(ipft) !------------------------------------------------------------------------------------! + return + end function dbh2bl + !=======================================================================================! + !=======================================================================================! + + + !=======================================================================================! + ! INVERSION OF DBH2BL ! + !=======================================================================================! + real function bl2dbh(bleaf,ipft) + use pft_coms, only : is_tropical & ! intent(in), lookup table + , rho & ! intent(in), lookup table + , dbh_crit & ! intent(in), lookup table + , hgt_max & ! intent(in), lookup table + , is_grass & ! intent(in) + , C2B & ! intent(in) + , b1Bl & ! intent(in), lookup table + , b2Bl ! ! intent(in), lookup table + use ed_misc_coms, only : igrass ! ! intent(in) + + implicit none + !----- Arguments --------------------------------------------------------------------! + real , intent(in) :: bleaf + integer, intent(in) :: ipft + !----- Local variables --------------------------------------------------------------! + real :: mdbh !------------------------------------------------------------------------------------! - ! Find maximum leaf biomass. ! - !------------------------------------------------------------------------------------! - dbh2bl = b1Bl(ipft) / C2B * mdbh ** b2Bl(ipft) + + mdbh = (bleaf * C2B / b1Bl(ipft) ) ** (1./b2Bl(ipft)) + + if (is_grass(ipft) .and. igrass==1) then + ! For grasses, limit maximum effective dbh by maximum height + bl2dbh = min(mdbh, h2dbh(hgt_max(ipft),ipft)) + else + bl2dbh = min(mdbh, dbh_crit(ipft)) + end if + + return + end function bl2dbh + !=======================================================================================! + !=======================================================================================! + + + + !=======================================================================================! + !=======================================================================================! + real function bl2h(bleaf,ipft) + use pft_coms, only: hgt_max ! ! intent(in), lookup table + + implicit none + !----- Arguments --------------------------------------------------------------------! + real , intent(in) :: bleaf + integer, intent(in) :: ipft !------------------------------------------------------------------------------------! + + !---Use existing allometric equations to convert leaves to height + bl2h = dbh2h(ipft,bl2dbh(bleaf,ipft)) + + if (bl2h > hgt_max(ipft)) then + bl2h = hgt_max(ipft) + end if return - end function dbh2bl + end function bl2h !=======================================================================================! !=======================================================================================! + !=======================================================================================! + !=======================================================================================! + real function h2bl(hite,ipft) + + implicit none + !----- Arguments --------------------------------------------------------------------! + real , intent(in) :: hite + integer, intent(in) :: ipft + !------------------------------------------------------------------------------------! + + !---Use existing allometric equations to convert height to leaves + h2bl = dbh2bl(h2dbh(hite,ipft), ipft) + + + return + end function h2bl + !=======================================================================================! + !=======================================================================================! @@ -205,15 +331,19 @@ end function dbh2bl !=======================================================================================! ! Canopy Area allometry from Dietze and Clark (2008). ! !---------------------------------------------------------------------------------------! - real function dbh2ca(dbh,sla,ipft) + real function dbh2ca(dbh,hite,sla,ipft) use ed_misc_coms, only : iallom ! ! intent(in) use pft_coms , only : dbh_crit & ! intent(in) + , hgt_max & ! intent(in) , is_tropical & ! intent(in) + , is_grass & ! intent(in) , b1Ca & ! intent(in) , b2Ca ! ! intent(in) + use ed_misc_coms, only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! real , intent(in) :: dbh + real , intent(in) :: hite real , intent(in) :: sla integer, intent(in) :: ipft !----- Internal variables -----------------------------------------------------------! @@ -224,7 +354,9 @@ real function dbh2ca(dbh,sla,ipft) loclai = 0.0 dbh2ca = 0.0 else - loclai = sla * dbh2bl(dbh,ipft) + + !!loclai = sla * dbh2bl(dbh,ipft) + loclai = sla * size2bl(dbh,hite,ipft) ! make this function generic to size, not just dbh select case (iallom) case (0) @@ -233,8 +365,11 @@ real function dbh2ca(dbh,sla,ipft) case default !----- Impose a maximum crown area. -------------------------------------------! - dbh2ca = b1Ca(ipft) * min(dbh,dbh_crit(ipft)) ** b2Ca(ipft) - + if (is_grass(ipft) .and. igrass==1) then + dbh2ca = b1Ca(ipft) * min(dbh,h2dbh(hgt_max(ipft),ipft) ) ** b2Ca(ipft) + else + dbh2ca = b1Ca(ipft) * min(dbh,dbh_crit(ipft) ) ** b2Ca(ipft) + end if end select end if @@ -280,7 +415,7 @@ end function dbh2vol !=======================================================================================! !=======================================================================================! - integer function dbh2krdepth(hgt,dbh,ipft,lsl) + integer function dbh2krdepth(hite,dbh,ipft,lsl) use ed_misc_coms, only : iallom ! ! intent(in) use grid_coms , only : nzg ! ! intent(in) use soil_coms , only : slz ! ! intent(in) @@ -288,7 +423,7 @@ integer function dbh2krdepth(hgt,dbh,ipft,lsl) , b2Rd ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! - real , intent(in) :: hgt + real , intent(in) :: hite real , intent(in) :: dbh integer, intent(in) :: ipft integer, intent(in) :: lsl @@ -305,7 +440,7 @@ integer function dbh2krdepth(hgt,dbh,ipft,lsl) !---------------------------------------------------------------------------------! ! Original ED-2.1 (I don't know the source for this equation, though). ! !---------------------------------------------------------------------------------! - volume = dbh2vol(hgt,dbh,ipft) + volume = dbh2vol(hite,dbh,ipft) root_depth = b1Rd(ipft) * volume ** b2Rd(ipft) case (1,2) @@ -313,7 +448,7 @@ integer function dbh2krdepth(hgt,dbh,ipft,lsl) ! This is just a test allometry, that imposes root depth to be 0.5 m for ! ! plants that are 0.15-m tall, and 5.0 m for plants that are 35-m tall. ! !---------------------------------------------------------------------------------! - root_depth = b1Rd(ipft) * hgt ** b2Rd(ipft) + root_depth = b1Rd(ipft) * hite ** b2Rd(ipft) end select @@ -369,56 +504,24 @@ end function h2crownbh - !=======================================================================================! - !=======================================================================================! - ! This subroutine finds the total above ground biomass corresponding to stems. ! - !---------------------------------------------------------------------------------------! - real function wood_biomass(bdead, bsapwood, pft) - use pft_coms, only: agf_bs ! ! intent(in) - - implicit none - !----- Arguments --------------------------------------------------------------------! - real , intent(in) :: bdead - real , intent(in) :: bsapwood - integer , intent(in) :: pft - !----- Local variables --------------------------------------------------------------! - real :: bstem - real :: absapwood - !------------------------------------------------------------------------------------! - - bstem = agf_bs(pft) * bdead - absapwood = agf_bs(pft) * bsapwood - wood_biomass = bstem + absapwood - return - end function wood_biomass - !=======================================================================================! - !=======================================================================================! - - - - - - !=======================================================================================! !=======================================================================================! ! This subroutine finds the total above ground biomass (wood + leaves) ! !---------------------------------------------------------------------------------------! - real function ed_biomass(bdead, balive, bleaf, pft, hite, bstorage, bsapwood) - + real function ed_biomass(bdead, bleaf, bsapwooda, ipft) + use pft_coms, only: agf_bs ! ! intent(in) + implicit none !----- Arguments --------------------------------------------------------------------! real , intent(in) :: bdead - real , intent(in) :: balive real , intent(in) :: bleaf - real , intent(in) :: hite - real , intent(in) :: bstorage - real , intent(in) :: bsapwood - integer , intent(in) :: pft + real , intent(in) :: bsapwooda + integer , intent(in) :: ipft !----- Local variables --------------------------------------------------------------! real :: bwood !------------------------------------------------------------------------------------! - bwood = wood_biomass(bdead, bsapwood, pft) + bwood = bsapwooda + (bdead * agf_bs(ipft)) ed_biomass = bleaf + bwood return @@ -444,8 +547,10 @@ end function ed_biomass ! didn't develop the allometry, but the original reference is in German...) ! !---------------------------------------------------------------------------------------! subroutine area_indices(nplant,bleaf,bdead,balive,dbh,hite,pft,sla,lai,wai,crown_area & - ,bsapwood) + ,bsapwooda) use pft_coms , only : is_tropical & ! intent(in) + , is_grass & ! intent(in) + , agf_bs & ! intent(in) , rho & ! intent(in) , C2B & ! intent(in) , dbh_crit & ! intent(in) @@ -454,6 +559,7 @@ subroutine area_indices(nplant,bleaf,bdead,balive,dbh,hite,pft,sla,lai,wai,crown use consts_coms , only : onethird & ! intent(in) , pi1 ! ! intent(in) use rk4_coms , only : ibranch_thermo ! ! intent(in) + use ed_misc_coms, only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! integer , intent(in) :: pft ! Plant functional type [ ---] @@ -461,7 +567,7 @@ subroutine area_indices(nplant,bleaf,bdead,balive,dbh,hite,pft,sla,lai,wai,crown real , intent(in) :: bleaf ! Specific leaf biomass [ kgC/plant] real , intent(in) :: bdead ! Specific structural [ kgC/plant] real , intent(in) :: balive ! Specific live tissue biomass [ kgC/plant] - real , intent(in) :: bsapwood ! Specific sapwood biomass [ kgC/plant] + real , intent(in) :: bsapwooda ! Specific sapwood biomass above grnd[ kgC/plant] real , intent(in) :: dbh ! Diameter at breast height [ cm] real , intent(in) :: hite ! Plant height [ m] real , intent(in) :: sla ! Specific leaf area [mēleaf/plant] @@ -482,9 +588,8 @@ subroutine area_indices(nplant,bleaf,bdead,balive,dbh,hite,pft,sla,lai,wai,crown !----- First, we compute the LAI ----------------------------------------------------! lai = bleaf * nplant * sla - !----- Find the crown area. ---------------------------------------------------------! - crown_area = min(1.0, nplant * dbh2ca(dbh,sla,pft)) + crown_area = min(1.0, nplant * dbh2ca(dbh,hite,sla,pft)) !------------------------------------------------------------------------------------! ! Here we check whether we need to compute the branch, stem, and effective ! diff --git a/ED/src/utils/ed_therm_lib.f90 b/ED/src/utils/ed_therm_lib.f90 index 028148afd..a8739b887 100644 --- a/ED/src/utils/ed_therm_lib.f90 +++ b/ED/src/utils/ed_therm_lib.f90 @@ -24,7 +24,7 @@ module ed_therm_lib ! to total biomass. The right hand side of the main equation accounts ! ! for the mass of insterstitial water and its ability to hold energy. ! ! + BDEAD - the structural wood biomass of the cohort in kgC/plant. ! - ! + BSAPWOOD - the sapwood biomass of the cohort, in kgC/plant. ! + ! + BSAPWOODA - the above ground sapwood biomass of the cohort, in kgC/plant. ! ! + NPLANTS - the number of plants per m2. ! ! + PFT - the plant functional type of the current cohort, which may serve ! ! for defining different parameterizations of specific heat capacity ! @@ -44,22 +44,23 @@ module ed_therm_lib ! energy storages on the land surface fluxes and radiative temperature. ! ! J. Geophys. Res., v. 112, doi: 10.1029/2006JD007425. ! !---------------------------------------------------------------------------------------! - subroutine calc_veg_hcap(bleaf,bdead,bsapwood,nplant,pft,leaf_hcap,wood_hcap) + subroutine calc_veg_hcap(bleaf,bdead,bsapwooda,nplant,pft,leaf_hcap,wood_hcap) use consts_coms , only : cliq ! ! intent(in) use pft_coms , only : c_grn_leaf_dry & ! intent(in) , wat_dry_ratio_grn & ! intent(in) , c_ngrn_biom_dry & ! intent(in) , wat_dry_ratio_ngrn & ! intent(in) , delta_c & ! intent(in) + , agf_bs & ! intent(in) , C2B & ! intent(in) , brf_wd ! ! intent(in) + use rk4_coms , only : ibranch_thermo ! ! intent(in) - use allometry , only : wood_biomass ! ! function implicit none !----- Arguments --------------------------------------------------------------------! real , intent(in) :: bleaf ! Biomass of leaves [kgC/plant] real , intent(in) :: bdead ! Biomass of structural wood [kgC/plant] - real , intent(in) :: bsapwood ! Biomass of sapwood [kgC/plant] + real , intent(in) :: bsapwooda ! Biomass of above ground sapwood[kgC/plant] real , intent(in) :: nplant ! Number of plants [ plant/m2] integer , intent(in) :: pft ! Plant functional type [ ----] real , intent(out) :: leaf_hcap ! Leaf heat capacity [ J/m2/K] @@ -82,10 +83,7 @@ subroutine calc_veg_hcap(bleaf,bdead,bsapwood,nplant,pft,leaf_hcap,wood_hcap) !----- Find branch/twig specific heat and biomass. -------------------------------! spheat_wood = (c_ngrn_biom_dry(pft) + wat_dry_ratio_ngrn(pft) * cliq) & / (1. + wat_dry_ratio_ngrn(pft)) + delta_c(pft) - - !----- Find the total branchwood biomass. ----------------------------------------! - bwood = brf_wd(pft) * wood_biomass(bdead, bsapwood,pft) - + bwood = brf_wd(pft) * (bsapwooda + bdead*agf_bs(pft)) end select !----- Find the leaf specific heat. -------------------------------------------------! diff --git a/ED/src/utils/fuse_fiss_utils.f90 b/ED/src/utils/fuse_fiss_utils.f90 index 19601a73e..2350cc98e 100644 --- a/ED/src/utils/fuse_fiss_utils.f90 +++ b/ED/src/utils/fuse_fiss_utils.f90 @@ -643,6 +643,7 @@ subroutine fuse_cohorts(csite,ipa, green_leaf_factor, lsl) , b1Ht & ! intent(in) , hgt_max & ! intent(in) , sla & ! intent(in) + , is_grass & ! intent(in) , hgt_ref ! ! intent(in) use fusion_fission_coms , only : fusetol_h & ! intent(in) , fusetol & ! intent(in) @@ -654,6 +655,7 @@ subroutine fuse_cohorts(csite,ipa, green_leaf_factor, lsl) use canopy_layer_coms , only : crown_mod ! ! intent(in) use allometry , only : dbh2h & ! function , dbh2bl ! ! function + use ed_misc_coms , only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! type(sitetype) , target :: csite ! Current site @@ -763,11 +765,19 @@ subroutine fuse_cohorts(csite,ipa, green_leaf_factor, lsl) ! leaves fully flushed, this is the same as adding the individual LAIs, ! ! but if they are not, we need to consider that LAI may grow... ! !------------------------------------------------------------------------! - lai_max = ( cpatch%nplant(recc) & - * dbh2bl(cpatch%dbh(recc),cpatch%pft(recc)) & - + cpatch%nplant(donc) & - * dbh2bl(cpatch%dbh(donc),cpatch%pft(donc))) & - * cpatch%sla(recc) + if (is_grass(cpatch%pft(donc)) .and. igrass==1) then + !--use actual bleaf for grass + lai_max = ( cpatch%nplant(recc) * cpatch%bleaf(recc) & + + cpatch%nplant(donc) * cpatch%bleaf(donc) ) & + * cpatch%sla(recc) + else + !--use dbh for trees + lai_max = ( cpatch%nplant(recc) & + * dbh2bl(cpatch%dbh(recc),cpatch%pft(recc)) & + + cpatch%nplant(donc) & + * dbh2bl(cpatch%dbh(donc),cpatch%pft(donc))) & + * cpatch%sla(recc) + end if !----- Checking the total size of this cohort before and after fusion. --! total_size = cpatch%nplant(donc) * ( cpatch%balive(donc) & @@ -777,7 +787,9 @@ subroutine fuse_cohorts(csite,ipa, green_leaf_factor, lsl) + cpatch%bdead(recc) & + cpatch%bstorage(recc) ) - + + + !------------------------------------------------------------------------! ! Five conditions must be met to allow two cohorts to be fused: ! ! 1. Both cohorts must have the same PFT; ! @@ -908,15 +920,19 @@ subroutine split_cohorts(cpatch, green_leaf_factor, lsl) use ed_state_vars , only : patchtype ! ! structure use pft_coms , only : q & ! intent(in), lookup table - , qsw ! ! intent(in), lookup table + , qsw & ! intent(in), lookup table + , is_grass ! ! intent(in) use fusion_fission_coms , only : lai_tol ! ! intent(in) use ed_max_dims , only : n_pft ! ! intent(in) use allometry , only : dbh2h & ! function , bd2dbh & ! function + , bl2dbh & ! function + , bl2h & ! function , dbh2bd ! ! function use ed_misc_coms , only : iqoutput & ! intent(in) , imoutput & ! intent(in) - , idoutput ! ! intent(in) + , idoutput & ! intent(in) + , igrass ! ! intent(in) use canopy_layer_coms , only : crown_mod ! ! intent(in) implicit none !----- Constants --------------------------------------------------------------------! @@ -1062,13 +1078,25 @@ subroutine split_cohorts(cpatch, green_leaf_factor, lsl) !---------------------------------------------------------------------------! !----- Tweaking bdead, to ensure carbon is conserved. ----------------------! - cpatch%bdead(ico) = cpatch%bdead(ico) * (1.-epsilon) - cpatch%dbh (ico) = bd2dbh(cpatch%pft(ico), cpatch%bdead(ico)) - cpatch%hite (ico) = dbh2h(cpatch%pft(ico), cpatch%dbh(ico)) - - cpatch%bdead(inew) = cpatch%bdead(inew) * (1.+epsilon) - cpatch%dbh (inew) = bd2dbh(cpatch%pft(inew), cpatch%bdead(inew)) - cpatch%hite (inew) = dbh2h(cpatch%pft(inew), cpatch%dbh(inew)) + if (is_grass(cpatch%pft(ico)) .and. igrass==1) then + !-- use bleaf for grass + cpatch%bleaf(ico) = cpatch%bleaf(ico) * (1.-epsilon) + cpatch%dbh (ico) = bl2dbh(cpatch%bleaf(ico), cpatch%pft(ico)) + cpatch%hite (ico) = bl2h(cpatch%bleaf(ico), cpatch%pft(ico)) + + cpatch%bleaf(inew) = cpatch%bleaf(inew) * (1.+epsilon) + cpatch%dbh (inew) = bl2dbh(cpatch%bleaf(inew), cpatch%pft(inew)) + cpatch%hite (inew) = bl2h(cpatch%bleaf(inew), cpatch%pft(inew)) + else + !-- use bdead for trees + cpatch%bdead(ico) = cpatch%bdead(ico) * (1.-epsilon) + cpatch%dbh (ico) = bd2dbh(cpatch%pft(ico), cpatch%bdead(ico)) + cpatch%hite (ico) = dbh2h(cpatch%pft(ico), cpatch%dbh(ico)) + + cpatch%bdead(inew) = cpatch%bdead(inew) * (1.+epsilon) + cpatch%dbh (inew) = bd2dbh(cpatch%pft(inew), cpatch%bdead(inew)) + cpatch%hite (inew) = dbh2h(cpatch%pft(inew), cpatch%dbh(inew)) + end if !---------------------------------------------------------------------------! end if @@ -1135,7 +1163,8 @@ subroutine clone_cohort(cpatch,isc,idt) cpatch%bdead(idt) = cpatch%bdead(isc) cpatch%bleaf(idt) = cpatch%bleaf(isc) cpatch%broot(idt) = cpatch%broot(isc) - cpatch%bsapwood(idt) = cpatch%bsapwood(isc) + cpatch%bsapwooda(idt) = cpatch%bsapwooda(isc) + cpatch%bsapwoodb(idt) = cpatch%bsapwoodb(isc) cpatch%phenology_status(idt) = cpatch%phenology_status(isc) cpatch%balive(idt) = cpatch%balive(isc) cpatch%lai(idt) = cpatch%lai(isc) @@ -1346,17 +1375,21 @@ end subroutine clone_cohort subroutine fuse_2_cohorts(cpatch,donc,recc, newn,green_leaf_factor, can_prss,lsl) use ed_state_vars , only : patchtype ! ! Structure use pft_coms , only : q & ! intent(in), lookup table - , qsw ! ! intent(in), lookup table + , qsw & ! intent(in), lookup table + , is_grass ! ! intent(in) use therm_lib , only : uextcm2tl & ! subroutine , qslif ! ! function use allometry , only : dbh2krdepth & ! function , bd2dbh & ! function + , bl2dbh & ! function + , bl2h & ! function , dbh2h ! ! function use ed_max_dims , only : n_mort ! ! intent(in) use ed_misc_coms , only : imoutput & ! intent(in) , iqoutput & ! intent(in) , idoutput & ! intent(in) - , ndcycle ! ! intent(in) + , ndcycle & ! intent(in) + , igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! type(patchtype) , target :: cpatch ! Current patch @@ -1394,12 +1427,23 @@ subroutine fuse_2_cohorts(cpatch,donc,recc, newn,green_leaf_factor, can_prss,lsl !----- Conserve carbon by calculating bdead first. ----------------------------------! - cpatch%bdead(recc) = ( cpatch%nplant(recc) * cpatch%bdead(recc) & - + cpatch%nplant(donc) * cpatch%bdead(donc) ) * newni - !----- Then get dbh and hite from bdead. --------------------------------------------! - cpatch%dbh(recc) = bd2dbh(cpatch%pft(recc), cpatch%bdead(recc)) - cpatch%hite(recc) = dbh2h(cpatch%pft(recc), cpatch%dbh(recc)) + if (is_grass(cpatch%pft(donc)) .and. igrass==1) then + !--use bleaf for grass + cpatch%bleaf(recc) = ( cpatch%nplant(recc) * cpatch%bleaf(recc) & + + cpatch%nplant(donc) * cpatch%bleaf(donc) ) * newni + cpatch%dbh(recc) = bl2dbh(cpatch%bleaf(recc), cpatch%pft(recc)) + cpatch%hite(recc) = bl2h (cpatch%bleaf(recc), cpatch%pft(recc)) + else + !--use bdead for trees + cpatch%bdead(recc) = ( cpatch%nplant(recc) * cpatch%bdead(recc) & + + cpatch%nplant(donc) * cpatch%bdead(donc) ) * newni + cpatch%dbh(recc) = bd2dbh(cpatch%pft(recc), cpatch%bdead(recc)) + cpatch%hite(recc) = dbh2h(cpatch%pft(recc), cpatch%dbh(recc)) + end if + + + !------------------------------------------------------------------------------------! @@ -1409,8 +1453,10 @@ subroutine fuse_2_cohorts(cpatch,donc,recc, newn,green_leaf_factor, can_prss,lsl + cpatch%nplant(donc) * cpatch%balive(donc) ) *newni cpatch%broot(recc) = ( cpatch%nplant(recc) * cpatch%broot(recc) & + cpatch%nplant(donc) * cpatch%broot(donc) ) *newni - cpatch%bsapwood(recc) = ( cpatch%nplant(recc) * cpatch%bsapwood(recc) & - + cpatch%nplant(donc) * cpatch%bsapwood(donc) ) *newni + cpatch%bsapwooda(recc) = ( cpatch%nplant(recc) * cpatch%bsapwooda(recc) & + + cpatch%nplant(donc) * cpatch%bsapwooda(donc) ) *newni + cpatch%bsapwoodb(recc) = ( cpatch%nplant(recc) * cpatch%bsapwoodb(recc) & + + cpatch%nplant(donc) * cpatch%bsapwoodb(donc) ) *newni cpatch%bstorage(recc) = ( cpatch%nplant(recc) * cpatch%bstorage(recc) & + cpatch%nplant(donc) * cpatch%bstorage(donc) ) * newni cpatch%bseeds(recc) = ( cpatch%nplant(recc) * cpatch%bseeds(recc) & @@ -4050,7 +4096,9 @@ subroutine patch_pft_size_profile(csite,ipa) , hgt_class ! ! intent(in) use allometry , only : dbh2bl ! ! intent(in) use ed_max_dims , only : n_pft ! ! intent(in) - use pft_coms , only : hgt_min ! ! intent(in) + use pft_coms , only : hgt_min & ! intent(in) + , is_grass ! ! intent(in) + use ed_misc_coms , only : igrass ! ! intent(in) implicit none !----- Arguments --------------------------------------------------------------------! type(sitetype) , target :: csite ! Current site @@ -4097,8 +4145,14 @@ subroutine patch_pft_size_profile(csite,ipa) !----- Find the potential (on-allometry) leaf area index. ------------------------! - lai_pot = cpatch%nplant(ico) * cpatch%sla(ico) & - * dbh2bl(cpatch%dbh(ico),ipft) + if (is_grass(ipft) .and. igrass==1) then + !--use actual bleaf for grass + lai_pot = cpatch%nplant(ico) * cpatch%sla(ico) * cpatch%bleaf(ico) + else + !--use dbh for trees + lai_pot = cpatch%nplant(ico) * cpatch%sla(ico) & + * dbh2bl(cpatch%dbh(ico),ipft) + end if !---------------------------------------------------------------------------------! diff --git a/ED/src/utils/update_derived_props.f90 b/ED/src/utils/update_derived_props.f90 index ce9ca5117..cac428415 100644 --- a/ED/src/utils/update_derived_props.f90 +++ b/ED/src/utils/update_derived_props.f90 @@ -116,10 +116,8 @@ subroutine update_patch_derived_props(csite,lsl,prss,ipa) !----- Compute the patch-level above-ground biomass csite%plant_ag_biomass(ipa) = csite%plant_ag_biomass(ipa) & - + ed_biomass(cpatch%bdead(ico),cpatch%balive(ico) & - ,cpatch%bleaf(ico),cpatch%pft(ico) & - ,cpatch%hite(ico),cpatch%bstorage(ico) & - ,cpatch%bsapwood(ico)) & + + ed_biomass(cpatch%bdead(ico),cpatch%bleaf(ico) & + ,cpatch%bsapwooda(ico),cpatch%pft(ico)) & * cpatch%nplant(ico) !------------------------------------------------------------------------------------! @@ -285,7 +283,6 @@ subroutine update_site_derived_props(cpoly,census_flag,isi) use ed_state_vars , only : polygontype & ! structure , sitetype & ! structure , patchtype ! ! structure - use allometry , only : ed_biomass ! ! function use consts_coms , only : pio4 ! ! intent(in) implicit none !----- Arguments -----------------------------------------------------------------------!