From 32dd172051dc3ace597db707bf30aba0610b7990 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 7 Oct 2016 16:52:56 -0700 Subject: [PATCH 01/11] Science fixes on sapwood and fine root respiration. First pass, no testing, compiles. --- .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 10 -- .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 114 ++++++++++-------- components/clm/src/ED/main/EDCLMLinkMod.F90 | 21 ---- components/clm/src/ED/main/EDTypesMod.F90 | 12 +- 4 files changed, 75 insertions(+), 82 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index fca32709d3..7c948f0b48 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -381,11 +381,6 @@ subroutine nan_cohort(cc_p) currentCohort%leaf_litter = nan ! leaf litter from phenology: KgC/m2 currentCohort%woody_turnover = nan ! amount of wood lost each day: kgC/indiv/year. Currently set to zero. - ! NITROGEN POOLS - currentCohort%livestemn = nan ! live stem nitrogen : KgN/invid - currentCohort%livecrootn = nan ! live coarse root nitrogen: KgN/invid - currentCohort%frootn = nan ! fine root nitrogen : KgN/invid - ! VARIABLES NEEDED FOR INTEGRATION currentCohort%dndt = nan ! time derivative of cohort size currentCohort%dhdt = nan ! time derivative of height @@ -1047,11 +1042,6 @@ subroutine copy_cohort( currentCohort,copyc ) n%livecroot_mr = o%livecroot_mr n%froot_mr = o%froot_mr - ! NITROGEN POOLS - n%livestemn = o%livestemn - n%livecrootn = o%livecrootn - n%frootn = o%frootn - ! ALLOCATION n%md = o%md n%leaf_md = o%leaf_md diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index ecad435267..f46d5e40ad 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -46,7 +46,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) use clm_varpar , only : nlevsoi, mxpft use clm_varctl , only : iulog use pftconMod , only : pftcon - use EDParamsMod , only : ED_val_grperc + use EDParamsMod , only : ED_val_grperc, & + ED_val_ag_biomass use EDSharedParamsMod , only : EDParamsShareInst use EDTypesMod , only : numpft_ed, dinc_ed use EDtypesMod , only : ed_patch_type, ed_cohort_type, ed_site_type, numpft_ed @@ -181,7 +182,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: an(cp_nclmax,mxpft,cp_nlevcan) ! net leaf photosynthesis (umol CO2/m**2/s) real(r8) :: an_av(cp_nclmax,mxpft,cp_nlevcan) ! net leaf photosynthesis (umol CO2/m**2/s) averaged over sun and shade leaves. real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) - real(r8) :: laican ! canopy sum of lai_z real(r8) :: vai ! leaf and steam area in ths layer. integer :: exitloop @@ -189,7 +189,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: tcsoi ! Temperature response function for root respiration. real(r8) :: tc ! Temperature response function for wood - real(r8) :: br ! Base rate of root respiration. (gC/gN/s) + real(r8) :: q10 ! temperature dependence of root respiration integer :: sunsha ! sun (1) or shaded (2) leaves... real(r8) :: dr(2) @@ -198,6 +198,26 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: gs_cohort real(r8) :: rscanopy real(r8) :: elai + + real(r8) :: live_agstem_n ! Live above-ground stem (sapwood) nitrogen content (gN/plant) + real(r8) :: live_bgstem_n ! Live below-ground stem (sapwood) nitrogen content (gN/plant) + real(r8) :: froot_n ! Fine root nitrogen content (gN/plant) + + ! Parameters + + ! ----------------------------------------------------------------------- + ! base maintenance respiration rate for plant tissues base_mr_20 + ! M. Ryan, 1991. Effects of climate change on plant respiration. + ! Ecological Applications, 1(2), 157-167. + ! Original expression is br = 0.0106 molC/(molN h) + ! Conversion by molecular weights of C and N gives 2.525e-6 gC/(gN s) + ! + ! Base rate is at 20C. Adjust to 25C using the CN Q10 = 1.5 + ! (gC/gN/s) + ! ------------------------------------------------------------------------ + + real(r8),parameter :: base_mr_20 = 2.525e-6_r8 + associate( & c3psn => pftcon%c3psn , & ! photosynthetic pathway: 0. = c4, 1. = c3 @@ -410,13 +430,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Leaf maintenance respiration to match the base rate used in CN ! but with the new temperature functions for C3 and C4 plants. ! - ! Base rate for maintenance respiration is from: - ! M. Ryan, 1991. Effects of climate change on plant respiration. - ! Ecological Applications, 1(2), 157-167. - ! Original expression is br = 0.0106 molC/(molN h) - ! Conversion by molecular weights of C and N gives 2.525e-6 gC/(gN s) - ! - ! Base rate is at 20C. Adjust to 25C using the CN Q10 = 1.5 ! ! CN respiration has units: g C / g N [leaf] / s. This needs to be ! converted from g C / g N [leaf] / s to umol CO2 / m**2 [leaf] / s @@ -853,60 +866,67 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) !------------------------------------------------------------------------------ ! Calculate Whole Plant Respiration (this doesn't really need to be in this iteration at all, surely?) ! Leaf respn needs to be in the sub-layer loop to account for changing N through canopy. - ! - ! base rate for maintenance respiration is from: - ! M. Ryan, 1991. Effects of climate change on plant respiration. - ! Ecological Applications, 1(2), 157-167. - ! Original expression is br = 0.0106 molC/(molN h) - ! Conversion by molecular weights of C and N gives 2.525e-6 gC/(gN s) !------------------------------------------------------------------------------ - br = 2.525e-6_r8 - leaf_frac = 1.0_r8/(currentCohort%canopy_trim + EDecophyscon%sapwood_ratio(currentCohort%pft) * & currentCohort%hite + pftcon%froot_leaf(currentCohort%pft)) - currentCohort%bsw = EDecophyscon%sapwood_ratio(currentCohort%pft) * currentCohort%hite * & - (currentCohort%balive + currentCohort%laimemory)*leaf_frac - currentCohort%livestemn = currentCohort%bsw / pftcon%leafcn(currentCohort%pft) - currentCohort%livestem_mr = 0._r8 - currentCohort%livecroot_mr = 0._r8 - if ( DEBUG ) write(iulog,*) 'EDPhoto 874 ', currentCohort%livecroot_mr - if ( DEBUG ) write(iulog,*) 'EDPhoto 875 ', currentCohort%livecrootn - - if (woody(FT) == 1) then + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! THIS CALCULATION SHOULD BE MOVED TO THE ALLOMETRY MODULE (RGK 10-8-2016) + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + currentCohort%bsw = EDecophyscon%sapwood_ratio(currentCohort%pft) * & + currentCohort%hite * (currentCohort%balive + currentCohort%laimemory)*leaf_frac + + ! Calculate the amount of nitrogen in the above and below ground + ! stem and root pools, used for maint resp + ! ------------------------------------------------------------------ + live_agstem_n = ED_val_ag_biomass * currentCohort%bsw / & + pftcon%leafcn(currentCohort%pft) + live_bgstem_n = (1.0_r8-ED_val_ag_biomass) * currentCohort%bsw / & + pftcon%leafcn(currentCohort%pft) + froot_n = currentCohort%br / leafcn(currentCohort%pft) + + ! Above ground Live stem MR + ! ------------------------------------------------------------------ + if (woody(ft) == 1) then tc = q10**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8) - currentCohort%livestem_mr = currentCohort%livestemn * br * tc !*currentPatch%btran_ft(currentCohort%pft) - currentCohort%livecroot_mr = currentCohort%livecrootn * br * tc !*currentPatch%btran_ft(currentCohort%pft) - + currentCohort%livestem_mr = live_agstem_n * base_mr_20 * tc !convert from gC /indiv/s-1 to kgC/indiv/s-1 - ! TODO: CHANGE THAT 1000 to 1000.0_r8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - currentCohort%livestem_mr = currentCohort%livestem_mr /1000 - currentCohort%livecroot_mr = currentCohort%livecroot_mr /1000 + currentCohort%livestem_mr = currentCohort%livestem_mr /1000.0_r8 else - tc = 1.0_r8 currentCohort%livestem_mr = 0._r8 - currentCohort%livecroot_mr = 0._r8 end if - if (pftcon%woody(currentCohort%pft) == 1) then - coarse_wood_frac = 0.5_r8 - else - coarse_wood_frac = 0.0_r8 - end if - ! Soil temperature. + ! Fine Root MR + ! ------------------------------------------------------------------ currentCohort%froot_mr = 0._r8 - do j = 1,nlevsoi tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) - !fine root respn. - currentCohort%froot_mr = currentCohort%froot_mr + (1.0_r8 - coarse_wood_frac) * & - currentCohort%br*br*tcsoi * currentPatch%rootfr_ft(ft,j)/leafcn(currentCohort%pft) - ! convert from gC/indiv/s-1 to kgC/indiv/s-1 - currentCohort%froot_mr = currentCohort%froot_mr /1000.0_r8 + currentCohort%froot_mr = currentCohort%froot_mr + & + froot_n * base_mr_20 * tcsoi * currentPatch%rootfr_ft(ft,j) enddo + ! convert from gC/indiv/s-1 to kgC/indiv/s-1 + currentCohort%froot_mr = currentCohort%froot_mr /1000.0_r8 + + + ! Coarse Root MR + ! ------------------------------------------------------------------ + if (woody(ft) == 1) then + currentCohort%livecroot_mr = 0._r8 + do j = 1,nlevsoi + ! Soil temperature used to adjust base rate of MR + tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) + currentCohort%livecroot_mr = currentCohort%livecroot_mr + & + live_bgstem_n * base_mr_20 * tcsoi * & + currentPatch%rootfr_ft(ft,j)* currentPatch%rootfr_ft(ft,j) + enddo + ! convert from gC/indiv/s-1 to kgC/indiv/s-1 + currentCohort%livecroot_mr = currentCohort%livecroot_mr /1000.0_r8 + else + currentCohort%livecroot_mr = 0._r8 + end if ! convert gpp and resp from umol/indiv/s-1 to kgC/indiv/s-1 = X * 12 *10-6 * 10-3 !currentCohort%resp_m = currentCohort%rd * 12.0E-9 diff --git a/components/clm/src/ED/main/EDCLMLinkMod.F90 b/components/clm/src/ED/main/EDCLMLinkMod.F90 index 1c7d6649f2..eab8da5aaf 100755 --- a/components/clm/src/ED/main/EDCLMLinkMod.F90 +++ b/components/clm/src/ED/main/EDCLMLinkMod.F90 @@ -175,27 +175,6 @@ subroutine ed_clm_link( this, bounds, nsites, sites, fcolumn, waterstate_inst, c do while(associated(currentCohort)) ft = currentCohort%pft - currentCohort%livestemn = currentCohort%bsw / pftcon%leafcn(currentCohort%pft) - - currentCohort%livecrootn = 0.0_r8 - - if (pftcon%woody(ft) == 1) then - coarse_wood_frac = 0.5_r8 - else - coarse_wood_frac = 0.0_r8 - end if - - if ( DEBUG ) then - write(iulog,*) 'EDCLMLink 618 ',currentCohort%livecrootn - write(iulog,*) 'EDCLMLink 619 ',currentCohort%br - write(iulog,*) 'EDCLMLink 620 ',coarse_wood_frac - write(iulog,*) 'EDCLMLink 621 ',pftcon%leafcn(ft) - endif - - currentCohort%livecrootn = currentCohort%br * coarse_wood_frac / pftcon%leafcn(ft) - - if ( DEBUG ) write(iulog,*) 'EDCLMLink 625 ',currentCohort%livecrootn - currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore currentCohort%treelai = tree_lai(currentCohort) diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 567c889c08..f9dff3af02 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -216,7 +216,9 @@ module EDTypesMod real(r8) :: resp_g ! Growth respiration: kgC/indiv/timestep real(r8) :: resp_m ! Maintenance respiration: kgC/indiv/timestep real(r8) :: livestem_mr ! Live stem maintenance respiration: kgC/indiv/s - real(r8) :: livecroot_mr ! Live coarse root maintenance respiration: kgC/indiv/s + ! (Above ground) + real(r8) :: livecroot_mr ! Live stem maintenance respiration: kgC/indiv/s + ! (below ground) real(r8) :: froot_mr ! Live fine root maintenance respiration: kgC/indiv/s ! ALLOCATION @@ -239,9 +241,11 @@ module EDTypesMod real(r8) :: fmort ! fire mortality n/year ! NITROGEN POOLS - real(r8) :: livestemn ! live stem nitrogen : KgN/invid - real(r8) :: livecrootn ! live coarse root nitrogen: KgN/invid - real(r8) :: frootn ! fine root nitrogen : KgN/invid + ! ---------------------------------------------------------------------------------- + ! Nitrogen pools are not prognostic in the current implementation. + ! They are diagnosed during photosynthesis using a simple C2N parameter. Local values + ! used in that routine. + ! ---------------------------------------------------------------------------------- ! GROWTH DERIVIATIVES real(r8) :: dndt ! time derivative of cohort size : n/year From 58e27deb4e33ed752d958b2effdfc81c8496b7ef Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 7 Oct 2016 17:02:50 -0700 Subject: [PATCH 02/11] fixed typo where MR of bg sapwood was double multiplied by bg root fraction --- components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index f46d5e40ad..79ed24f3a0 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -920,7 +920,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%livecroot_mr = currentCohort%livecroot_mr + & live_bgstem_n * base_mr_20 * tcsoi * & - currentPatch%rootfr_ft(ft,j)* currentPatch%rootfr_ft(ft,j) + currentPatch%rootfr_ft(ft,j) enddo ! convert from gC/indiv/s-1 to kgC/indiv/s-1 currentCohort%livecroot_mr = currentCohort%livecroot_mr /1000.0_r8 From 4416814360f4f10679d76a026fc3a2a8f02cf0e0 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 10 Oct 2016 16:42:22 -0700 Subject: [PATCH 03/11] Added several AR diagnostics. Also changed the units of dark-respiration to be consistent with the other maintenance pools. However, the unit of kgC/plant/sec is likely to be a very very small number that is wasting precision. Contemplating changing units to umol/sec or kgC/yr. --- .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 6 +- .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 44 +++++---- components/clm/src/ED/main/EDTypesMod.F90 | 2 +- components/clm/src/ED/main/HistoryIOMod.F90 | 96 ++++++++++++++++++- 4 files changed, 123 insertions(+), 25 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 7c948f0b48..d35afdb1d8 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -361,7 +361,7 @@ subroutine nan_cohort(cc_p) !RESPIRATION - currentCohort%rd = nan + currentCohort%rdark = nan currentCohort%resp_m = nan ! Maintenance respiration. kGC/cohort/year currentCohort%resp_g = nan ! Growth respiration. kGC/cohort/year currentCohort%livestem_mr = nan ! Live stem maintenance respiration. kgC/indiv/s-1 @@ -418,7 +418,7 @@ subroutine zero_cohort(cc_p) currentCohort%NV = 0 currentCohort%status_coh = 0 - currentCohort%rd = 0._r8 + currentCohort%rdark = 0._r8 currentCohort%resp_m = 0._r8 currentCohort%resp_g = 0._r8 currentCohort%livestem_mr = 0._r8 @@ -1035,7 +1035,7 @@ subroutine copy_cohort( currentCohort,copyc ) n%npp_store = o%npp_store !RESPIRATION - n%rd = o%rd + n%rdark = o%rdark n%resp_m = o%resp_m n%resp_g = o%resp_g n%livestem_mr = o%livestem_mr diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 79ed24f3a0..35334fa7f5 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -205,6 +205,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Parameters + ! Conversion factor umols of Carbon -> kg of Carbon (1 mol = 12g) + real(r8), parameter :: umolC_to_kgC = 12.0E-9_r8 + + ! Conversion factor: grams per kilograms + real(r8), parameter :: g_per_kg = 1000.0_r8 + ! ----------------------------------------------------------------------- ! base maintenance respiration rate for plant tissues base_mr_20 ! M. Ryan, 1991. Effects of climate change on plant respiration. @@ -437,7 +443,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Then scale this value at the top of the canopy for canopy depth lmr25top(FT) = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) - lmr25top(FT) = lmr25top(FT) * lnc(FT) / 12.e-06_r8 + lmr25top(FT) = lmr25top(FT) * lnc(FT) / (umolC_to_kgC * g_per_kg) end do !FT @@ -801,7 +807,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%npp_tstep = 0.0_r8 currentCohort%resp_tstep = 0.0_r8 currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rd = 0.0_r8 + currentCohort%rdark = 0.0_r8 currentCohort%resp_m = 0.0_r8 ! Select canopy layer and PFT. @@ -824,16 +830,16 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%gpp_tstep = sum(currentPatch%psn_z(cl,ft,1:currentCohort%nv-1) * & currentPatch%elai_profile(cl,ft,1:currentCohort%nv-1)) * tree_area - currentCohort%rd = sum(lmr_z(cl,ft,1:currentCohort%nv-1) * & + currentCohort%rdark = sum(lmr_z(cl,ft,1:currentCohort%nv-1) * & currentPatch%elai_profile(cl,ft,1:currentCohort%nv-1)) * tree_area currentCohort%gscan = sum((1.0_r8/(rs_z(cl,ft,1:currentCohort%nv-1)+bc_in(s)%rb_pa(ifp)))) * tree_area - currentCohort%ts_net_uptake(1:currentCohort%nv) = an_av(cl,ft,1:currentCohort%nv) * 12E-9 * dtime + currentCohort%ts_net_uptake(1:currentCohort%nv) = an_av(cl,ft,1:currentCohort%nv) * umolC_to_kgC * dtime else currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rd = 0.0_r8 + currentCohort%rdark = 0.0_r8 currentCohort%gscan = 0.0_r8 currentCohort%ts_net_uptake(:) = 0.0_r8 @@ -858,15 +864,14 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%gpp_tstep = currentCohort%gpp_tstep + currentPatch%psn_z(cl,ft,currentCohort%nv) * & currentPatch%elai_profile(cl,ft,currentCohort%nv) * laifrac * tree_area - if ( DEBUG ) write(iulog,*) 'EDPhoto 843 ', currentCohort%rd + if ( DEBUG ) write(iulog,*) 'EDPhoto 843 ', currentCohort%rdark - currentCohort%rd = currentCohort%rd + lmr_z(cl,ft,currentCohort%nv) * & + + currentCohort%rdark = currentCohort%rdark + lmr_z(cl,ft,currentCohort%nv) * & currentPatch%elai_profile(cl,ft,currentCohort%nv) * laifrac * tree_area - !------------------------------------------------------------------------------ - ! Calculate Whole Plant Respiration (this doesn't really need to be in this iteration at all, surely?) - ! Leaf respn needs to be in the sub-layer loop to account for changing N through canopy. - !------------------------------------------------------------------------------ + ! Convert dark respiration from umol/plant/s to kgC/plant/s + currentCohort%rdark = currentCohort%rdark * umolC_to_kgC leaf_frac = 1.0_r8/(currentCohort%canopy_trim + EDecophyscon%sapwood_ratio(currentCohort%pft) * & currentCohort%hite + pftcon%froot_leaf(currentCohort%pft)) @@ -887,6 +892,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) pftcon%leafcn(currentCohort%pft) froot_n = currentCohort%br / leafcn(currentCohort%pft) + + !------------------------------------------------------------------------------ + ! Calculate Whole Plant Respiration (this doesn't really need to be in this iteration at all, surely?) + ! Leaf respn needs to be in the sub-layer loop to account for changing N through canopy. + !------------------------------------------------------------------------------ + ! Above ground Live stem MR ! ------------------------------------------------------------------ if (woody(ft) == 1) then @@ -928,20 +939,19 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%livecroot_mr = 0._r8 end if - ! convert gpp and resp from umol/indiv/s-1 to kgC/indiv/s-1 = X * 12 *10-6 * 10-3 - !currentCohort%resp_m = currentCohort%rd * 12.0E-9 + ! convert gpp from umol/indiv/s-1 to kgC/indiv/s-1 = X * 12 *10-6 * 10-3 if ( DEBUG ) write(iulog,*) 'EDPhoto 904 ', currentCohort%resp_m - if ( DEBUG ) write(iulog,*) 'EDPhoto 905 ', currentCohort%rd + if ( DEBUG ) write(iulog,*) 'EDPhoto 905 ', currentCohort%rdark if ( DEBUG ) write(iulog,*) 'EDPhoto 906 ', currentCohort%livestem_mr if ( DEBUG ) write(iulog,*) 'EDPhoto 907 ', currentCohort%livecroot_mr if ( DEBUG ) write(iulog,*) 'EDPhoto 908 ', currentCohort%froot_mr - currentCohort%gpp_tstep = currentCohort%gpp_tstep * 12.0E-9 + currentCohort%gpp_tstep = currentCohort%gpp_tstep * umolC_to_kgC ! add on whole plant respiration values in kgC/indiv/s-1 currentCohort%resp_m = currentCohort%livestem_mr + currentCohort%livecroot_mr + currentCohort%froot_mr ! no drought response * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft)*pftcon%resp_drought_response(FT)) - currentCohort%resp_m = currentCohort%resp_m + currentCohort%rd * 12.0E-9 !this was already corrected fo BTRAN + currentCohort%resp_m = currentCohort%resp_m + currentCohort%rdark ! convert from kgC/indiv/s to kgC/indiv/timestep currentCohort%resp_m = currentCohort%resp_m * dtime @@ -951,7 +961,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) if ( DEBUG ) write(iulog,*) 'EDPhoto 912 ', currentCohort%resp_tstep if ( DEBUG ) write(iulog,*) 'EDPhoto 913 ', currentCohort%resp_m - currentCohort%resp_g = ED_val_grperc(1) * (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) + currentCohort%resp_g = ED_val_grperc(1) * (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) currentCohort%resp_tstep = currentCohort%resp_m + currentCohort%resp_g ! kgC/indiv/ts currentCohort%npp_tstep = currentCohort%gpp_tstep - currentCohort%resp_tstep ! kgC/indiv/ts diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index f9dff3af02..4af5b32b57 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -212,7 +212,7 @@ module EDTypesMod real(r8) :: year_net_uptake(cp_nlevcan) ! Net uptake of leaf layers: kgC/m2/year ! RESPIRATION COMPONENTS - real(r8) :: rd ! Dark respiration: umol/indiv/s + real(r8) :: rdark ! Dark respiration: kgC/indiv/s real(r8) :: resp_g ! Growth respiration: kgC/indiv/timestep real(r8) :: resp_m ! Maintenance respiration: kgC/indiv/timestep real(r8) :: livestem_mr ! Live stem maintenance respiration: kgC/indiv/s diff --git a/components/clm/src/ED/main/HistoryIOMod.F90 b/components/clm/src/ED/main/HistoryIOMod.F90 index 6faa0d9eac..c746c33d30 100644 --- a/components/clm/src/ED/main/HistoryIOMod.F90 +++ b/components/clm/src/ED/main/HistoryIOMod.F90 @@ -119,6 +119,13 @@ Module HistoryIOMod integer, private :: ih_m4_si_scpf integer, private :: ih_m5_si_scpf + integer, private :: ih_ar_si_scpf + integer, private :: ih_ar_grow_si_scpf + integer, private :: ih_ar_maint_si_scpf + integer, private :: ih_ar_darkm_si_scpf + integer, private :: ih_ar_agsapm_si_scpf + integer, private :: ih_ar_crootm_si_scpf + integer, private :: ih_ar_frootm_si_scpf ! The number of variable dim/kind types we have defined (static) integer, parameter :: n_iovar_dk = 6 @@ -627,7 +634,9 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) use EDtypesMod , only : ed_site_type, & ed_cohort_type, & ed_patch_type, & - AREA + AREA, & + sclass_ed, & + nlevsclass_ed ! Arguments class(fates_hio_interface_type) :: this integer , intent(in) :: nc ! clump index @@ -644,7 +653,9 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) integer :: io_soipa integer :: lb1,ub1,lb2,ub2 ! IO array bounds for the calling thread integer :: ivar ! index of IO variable object vector - + integer :: ft ! functional type index + integer :: scpf ! index of the size-class x pft bin + integer :: sc ! index of the size-class bin real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column @@ -660,7 +671,15 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) hio_aresp_pa => this%hvars(ih_aresp_pa)%r81d, & hio_maint_resp_pa => this%hvars(ih_maint_resp_pa)%r81d, & hio_growth_resp_pa => this%hvars(ih_growth_resp_pa)%r81d, & - hio_npp_si => this%hvars(ih_npp_si)%r81d ) + hio_npp_si => this%hvars(ih_npp_si)%r81d, & + hio_ar_si_scpf => this%hvars(ih_ar_si_scpf)%r82d, & + hio_ar_grow_si_scpf => this%hvars(ih_ar_grow_si_scpf)%r82d, & + hio_ar_maint_si_scpf => this%hvars(ih_ar_maint_si_scpf)%r82d, & + hio_ar_agsapm_si_scpf => this%hvars(ih_ar_agsapm_si_scpf)%r82d, & + hio_ar_darkm_si_scpf => this%hvars(ih_ar_darkm_si_scpf)%r82d, & + hio_ar_crootm_si_scpf => this%hvars(ih_ar_crootm_si_scpf)%r82d, & + hio_ar_frootm_si_scpf => this%hvars(ih_ar_frootm_si_scpf)%r82d ) + ! Flush the relevant history variables call this%flush_hvars(nc,upfreq_in=2) @@ -706,6 +725,39 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! map ed cohort-level npp fluxes to clm column fluxes hio_npp_si(io_si) = hio_npp_si(io_si) + ccohort%npp_tstep * n_perm2 * 1.e3_r8 /dt_tstep + ! Calculate index for the scpf class + sc = count(ccohort%dbh-sclass_ed.ge.0.0) + scpf = (ft-1)*nlevsclass_ed+sc + + ! Total AR (kgC/m2/yr) = (kgC/plant/step) / (s/step) * (plant/m2) * (s/yr) + hio_ar_si_scpf(io_si,scpf) = hio_ar_si_scpf(io_si,scpf) + & + (ccohort%resp_tstep/dt_tstep) * n_perm2 * daysecs * yeardays + + ! Growth AR (kgC/m2/yr) + hio_ar_grow_si_scpf(io_si,scpf) = hio_ar_grow_si_scpf(io_si,scpf) + & + (ccohort%resp_g/dt_tstep) * n_perm2 * daysecs * yeardays + + ! Maint AR (kgC/m2/yr) + hio_ar_maint_si_scpf(io_si,scpf) = hio_ar_maint_si_scpf(io_si,scpf) + & + (ccohort%resp_m/dt_tstep) * n_perm2 * daysecs * yeardays + + ! Maintenance AR partition variables are stored as rates (kgC/plant/s) + ! (kgC/m2/yr) = (kgC/plant/s) * (plant/m2) * (s/yr) + hio_ar_agsapm_si_scpf(io_si,scpf) = hio_ar_agsapm_si_scpf(io_si,scpf) + & + ccohort%livestem_mr * n_perm2 * daysecs * yeardays + + ! (kgC/m2/yr) = (kgC/plant/s) * (plant/m2) * (s/yr) + hio_ar_darkm_si_scpf(io_si,scpf) = hio_ar_darkm_si_scpf(io_si,scpf) + & + ccohort%rdark * n_perm2 * daysecs * yeardays + + ! (kgC/m2/yr) = (kgC/plant/s) * (plant/m2) * (s/yr) + hio_ar_crootm_si_scpf(io_si,scpf) = hio_ar_crootm_si_scpf(io_si,scpf) + & + ccohort%livecroot_mr * n_perm2 * daysecs * yeardays + + ! (kgC/m2/yr) = (kgC/plant/s) * (plant/m2) * (s/yr) + hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + & + ccohort%froot_mr * n_perm2 * daysecs * yeardays + endif ccohort => ccohort%taller @@ -1028,7 +1080,6 @@ subroutine define_history_vars(this,callstep,nvar) avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & upfreq=1, ivar=ivar,callstep=callstep, index = ih_npp_totl_si_scpf ) - call this%set_history_var(vname='NPP_LEAF_SCPF',units='kgC/m2/yr', & long='NPP flux into leaves', use_default='inactive', & avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & @@ -1109,6 +1160,43 @@ subroutine define_history_vars(this,callstep,nvar) avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & upfreq=1, ivar=ivar,callstep=callstep, index = ih_m5_si_scpf ) + ! Size structured diagnostics that require rapid updates (upfreq=2) + + call this%set_history_var(vname='AR_SCPF',units = 'kgC/m2/yr', & + long='total autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_si_scpf ) + + call this%set_history_var(vname='AR_GROW_SCPF',units = 'kgC/m2/yr', & + long='growth autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_grow_si_scpf ) + + call this%set_history_var(vname='AR_MAINT_SCPF',units = 'kgC/m2/yr', & + long='maintenance autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_maint_si_scpf ) + + call this%set_history_var(vname='AR_DARKM_SCPF',units = 'kgC/m2/yr', & + long='dark portion of maintenance autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_darkm_si_scpf ) + + call this%set_history_var(vname='AR_AGSAPM_SCPF',units = 'kgC/m2/yr', & + long='above-ground sapwood maintenance autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_agsapm_si_scpf ) + + call this%set_history_var(vname='AR_CROOTM_SCPF',units = 'kgC/m2/yr', & + long='below-ground sapwood maintenance autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_crootm_si_scpf ) + + call this%set_history_var(vname='AR_FROOTM_SCPF',units = 'kgC/m2/yr', & + long='fine root maintenance autotrophic respiration per m2 per year',use_default='inactive',& + avgflag='A', vtype='SI_SCPF_R8',hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=2, ivar=ivar,callstep=callstep, index = ih_ar_frootm_si_scpf ) + ! CARBON BALANCE VARIABLES THAT DEPEND ON HLM BGC INPUTS From 0ca158e8c36aaefca58c4b3738e4a259f47b5200 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 11 Oct 2016 21:22:20 -0700 Subject: [PATCH 04/11] minor fix to update the pft index during the scpf output for photosynthesis --- components/clm/src/ED/main/HistoryIOMod.F90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/components/clm/src/ED/main/HistoryIOMod.F90 b/components/clm/src/ED/main/HistoryIOMod.F90 index c746c33d30..85c51ffd74 100644 --- a/components/clm/src/ED/main/HistoryIOMod.F90 +++ b/components/clm/src/ED/main/HistoryIOMod.F90 @@ -710,6 +710,11 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) if ( .not. ccohort%isnew ) then + ! Calculate index for the scpf class + ft = ccohort%pft + sc = count(ccohort%dbh-sclass_ed.ge.0.0) + scpf = (ft-1)*nlevsclass_ed+sc + ! scale up cohort fluxes to their patches hio_npp_pa(io_pa) = hio_npp_pa(io_pa) + & ccohort%npp_tstep * 1.e3_r8 * n_density / dt_tstep @@ -725,9 +730,6 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! map ed cohort-level npp fluxes to clm column fluxes hio_npp_si(io_si) = hio_npp_si(io_si) + ccohort%npp_tstep * n_perm2 * 1.e3_r8 /dt_tstep - ! Calculate index for the scpf class - sc = count(ccohort%dbh-sclass_ed.ge.0.0) - scpf = (ft-1)*nlevsclass_ed+sc ! Total AR (kgC/m2/yr) = (kgC/plant/step) / (s/step) * (plant/m2) * (s/yr) hio_ar_si_scpf(io_si,scpf) = hio_ar_si_scpf(io_si,scpf) + & From 4ecda622d404a3099222761294c3e2ad4f6d2a43 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 11 Oct 2016 21:37:58 -0700 Subject: [PATCH 05/11] minor fix, livecrootn was still in a print statement that I overlooked. I removed it. --- components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 index 5e3110973a..5b8ad923cf 100755 --- a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 @@ -711,8 +711,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ft = currentCohort%pft - if ( DEBUG ) write(fates_log(),*) 'canopy_summarization 732 ',currentCohort%livecrootn - currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore currentCohort%treelai = tree_lai(currentCohort) From d0150fe7184d602d480a3c95f0fd6c09292f8116 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 13 Oct 2016 14:51:34 -0700 Subject: [PATCH 06/11] fixed hard-coded pft index in grperc found in photosynthesis --- components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 | 2 +- components/clm/src/utils/clmfates_interfaceMod.F90 | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 18dd7e6d77..9671570320 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -960,7 +960,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) if ( DEBUG ) write(iulog,*) 'EDPhoto 912 ', currentCohort%resp_tstep if ( DEBUG ) write(iulog,*) 'EDPhoto 913 ', currentCohort%resp_m - currentCohort%resp_g = ED_val_grperc(1) * (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) + currentCohort%resp_g = ED_val_grperc(ft) * (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) currentCohort%resp_tstep = currentCohort%resp_m + currentCohort%resp_g ! kgC/indiv/ts currentCohort%npp_tstep = currentCohort%gpp_tstep - currentCohort%resp_tstep ! kgC/indiv/ts diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 86fd070b4b..0fbc7dd51c 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -1063,7 +1063,6 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, & use perf_mod , only : t_startf, t_stopf use PatchType , only : patch use quadraticMod , only : quadratic - use EDParamsMod , only : ED_val_grperc use EDSharedParamsMod , only : EDParamsShareInst use EDTypesMod , only : numpft_ed, dinc_ed use EDtypesMod , only : ed_patch_type, ed_cohort_type, ed_site_type, numpft_ed, numPatchesPerCol From 15824a24f34ec9c937577d66aa5aaabeb145b17b Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 13 Oct 2016 18:12:53 -0700 Subject: [PATCH 07/11] added some more global constants and uses in photosynthesis. Also, fixed a unit conversion error in maintenance respiration. --- .../ED/biogeochem/EDCanopyStructureMod.F90 | 1 - .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 212 ++++++++++-------- .../clm/src/ED/main/FatesConstantsMod.F90 | 29 +++ 3 files changed, 142 insertions(+), 100 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 index 5b8ad923cf..8ad6b71c79 100755 --- a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 @@ -675,7 +675,6 @@ subroutine canopy_summarization( nsites, sites, bc_in ) integer :: ft ! plant functional type integer :: ifp integer :: patchn ! identification number for each patch. - real(r8) :: coarse_wood_frac real(r8) :: canopy_leaf_area ! total amount of leaf area in the vegetated area. m2. !---------------------------------------------------------------------- diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 9671570320..79b4897728 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -9,14 +9,14 @@ module EDPhotosynthesisMod ! ! !USES: ! - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_log_mod , only : errMsg => shr_log_errMsg - use abortutils , only : endrun - use clm_varctl , only : iulog + + use abortutils, only : endrun + use FatesGlobals, only : fates_log + use FatesConstantsMod, only : r8 => fates_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg + implicit none private - ! - ! PUBLIC MEMBER FUNCTIONS: @@ -39,22 +39,39 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! a multi-layer canopy ! ! !USES: - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun - use clm_varcon , only : rgas, tfrz, namep - use clm_varpar , only : nlevsoi, mxpft - use clm_varctl , only : iulog - use pftconMod , only : pftcon + use clm_varpar , only : mxpft ! THIS WILL BE DEPRECATED WHEN PARAMETER + ! READS ARE REFACTORED (RGK 10-13-2016) + use pftconMod , only : pftcon ! THIS WILL BE DEPRECATED WHEN PARAMETER + ! READS ARE REFACTORED (RGK 10-13-2016) use EDParamsMod , only : ED_val_grperc, & ED_val_ag_biomass use EDSharedParamsMod , only : EDParamsShareInst - use EDTypesMod , only : numpft_ed, dinc_ed - use EDtypesMod , only : ed_patch_type, ed_cohort_type, ed_site_type, numpft_ed + use EDTypesMod , only : numpft_ed, & + dinc_ed, & + ed_patch_type, & + ed_cohort_type, & + ed_site_type, & + numpft_ed, & + numpatchespercol, & + cp_numlevsoil, & + cp_nlevcan, & + cp_nclmax + use EDEcophysContype , only : EDecophyscon - use FatesInterfaceMod , only : bc_in_type,bc_out_type - use EDtypesMod , only : numpatchespercol, cp_nlevcan, cp_nclmax + use FatesInterfaceMod , only : bc_in_type, & + bc_out_type + use EDCanopyStructureMod,only: calc_areaindex + + use FatesConstantsMod, only : umolC_to_kgC, & ! micromole conversion to kgC + g_per_kg, & ! number of grams per kg + mg_per_g, & ! number of miligrams per g + sec_per_min, & ! seconds per minute (60!) + rgas, & ! universal gas constant + + tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atmosphere ! ! !ARGUMENTS: @@ -149,7 +166,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) integer :: c,CL,f,s,iv,j,ps,ft,ifp ! indices integer :: NCL_p ! number of canopy layers in patch real(r8) :: cf ! s m**2/umol -> s/m - real(r8) :: rsmax0 ! maximum stomatal resistance [s/m] + real(r8) :: gb ! leaf boundary layer conductance (m/s) real(r8) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) real(r8) :: cs ! CO2 partial pressure at leaf surface (Pa) @@ -191,7 +208,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: q10 ! temperature dependence of root respiration integer :: sunsha ! sun (1) or shaded (2) leaves... - real(r8) :: dr(2) real(r8) :: coarse_wood_frac ! amount of woody biomass that is coarse... real(r8) :: tree_area real(r8) :: gs_cohort @@ -203,15 +219,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: froot_n ! Fine root nitrogen content (gN/plant) ! Parameters - - ! Conversion factor umols of Carbon -> kg of Carbon (1 mol = 12g) - real(r8), parameter :: umolC_to_kgC = 12.0E-9_r8 - - ! Conversion factor: grams per kilograms - real(r8), parameter :: g_per_kg = 1000.0_r8 - ! ----------------------------------------------------------------------- - ! base maintenance respiration rate for plant tissues base_mr_20 + ! Base maintenance respiration rate for plant tissues base_mr_20 ! M. Ryan, 1991. Effects of climate change on plant respiration. ! Ecological Applications, 1(2), 157-167. ! Original expression is br = 0.0106 molC/(molN h) @@ -222,7 +231,15 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! ------------------------------------------------------------------------ real(r8),parameter :: base_mr_20 = 2.525e-6_r8 - + + ! maximum stomatal resistance [s/m] + real(r8),parameter :: rsmax0 = 2.e4_r8 + + ! First guess on ratio between intracellular co2 and the atmosphere + ! an iterator converges on actual + real(r8),parameter :: init_a2l_co2_c3 = 0.7 + real(r8),parameter :: init_a2l_co2_c4 = 0.4 + associate( & c3psn => pftcon%c3psn , & ! photosynthetic pathway: 0. = c4, 1. = c3 @@ -231,11 +248,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) woody => pftcon%woody , & ! Is vegetation woody or not? fnitr => pftcon%fnitr , & ! foliage nitrogen limitation factor (-) leafcn => pftcon%leafcn , & ! leaf C:N (gC/gN) + frootcn => pftcon%frootcn , & ! froot C:N (gc/gN) bb_slope => EDecophyscon%BB_slope ) ! slope of BB relationship - ! Assign local pointers to derived type members (gridcell-level) - dr(1) = 0.025_r8; dr(2) = 0.015_r8 - ! Peter Thornton: 3/13/09 ! Q10 was originally set to 2.0, an arbitrary choice, but reduced to 1.5 as part of the tuning ! to improve seasonal cycle of atmospheric CO2 concentration in global @@ -252,7 +267,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) act25 = 3.6_r8 !umol/mgRubisco/min ! Convert rubisco activity units from umol/mgRubisco/min -> umol/gRubisco/s - act25 = act25 * 1000.0_r8 / 60.0_r8 + act25 = act25 * mg_per_g / sec_per_min ! Activation energy, from: ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 @@ -399,12 +414,11 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 650 (RGK) if (nint(c3psn(FT)) == 1)then - ci(:,FT,:) = 0.7_r8 * bc_in(s)%cair_pa(ifp) + ci(:,FT,:) = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) else - ci(:,FT,:) = 0.4_r8 * bc_in(s)%cair_pa(ifp) + ci(:,FT,:) = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) end if - ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) lnc(FT) = 1._r8 / (slatop(FT) * leafcn(FT)) @@ -416,7 +430,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Here use a factor "1.67", from Medlyn et al (2002) Plant, Cell and Environment 25:1167-1179 !RF - copied this from the CLM trunk code, but where did it come from, and how can we make these consistant? - !jmax25top(FT) = (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrz),11._r8),35._r8)) * vcmax25top(FT) + !jmax25top(FT) = (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(FT) jmax25top(FT) = 1.67_r8 * vcmax25top(FT) tpu25top(FT) = 0.167_r8 * vcmax25top(FT) @@ -454,7 +468,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) do iv = 1, currentPatch%nrad(CL,FT) if(currentPatch%canopy_area_profile(CL,FT,iv)>0._r8.and.currentPatch%present(CL,FT) /= 1)then - write(iulog,*) 'CF: issue with present structure',CL,FT,iv, & + write(fates_log(),*) 'CF: issue with present structure',CL,FT,iv, & currentPatch%canopy_area_profile(CL,FT,iv),currentPatch%present(CL,FT), & currentPatch%nrad(CL,FT),currentPatch%ncl_p,cp_nclmax currentPatch%present(CL,FT) = 1 @@ -539,10 +553,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Leaf-level photosynthesis and stomatal conductance !==============================================================================! - rsmax0 = 2.e4_r8 - ! Leaf boundary layer conductance, umol/m**2/s + ! THESE HARD CODED CONVERSIONS NEED TO BE CALLED FROM GLOBAL CONSTANTS (RGK 10-13-2016) cf = bc_in(s)%forc_pbot/(rgas*1.e-3_r8*bc_in(s)%tgcm_pa(ifp))*1.e06_r8 gb = 1._r8/bc_in(s)%rb_pa(ifp) gb_mol = gb * cf @@ -561,7 +574,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if if(currentPatch%present(CL,FT) == 1)then ! are there any leaves of this pft in this layer? do iv = 1, currentPatch%nrad(CL,FT) - if ( DEBUG ) write(iulog,*) 'EDphoto 581 ',currentPatch%ed_parsun_z(CL,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDphoto 581 ',currentPatch%ed_parsun_z(CL,ft,iv) if (currentPatch%ed_parsun_z(CL,FT,iv) <= 0._r8) then ! night time ac = 0._r8 @@ -576,12 +589,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) else ! day time !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) - if ( DEBUG ) write(iulog,*) 'EDphot 594 ',currentPatch%ed_laisun_z(CL,ft,iv) - if ( DEBUG ) write(iulog,*) 'EDphot 595 ',currentPatch%ed_laisha_z(CL,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDphot 594 ',currentPatch%ed_laisun_z(CL,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDphot 595 ',currentPatch%ed_laisha_z(CL,ft,iv) if(currentPatch%ed_laisun_z(CL,ft,iv)+currentPatch%ed_laisha_z(cl,ft,iv) > 0._r8)then - if ( DEBUG ) write(iulog,*) '600 in laisun, laisha loop ' + if ( DEBUG ) write(fates_log(),*) '600 in laisun, laisha loop ' !Loop aroun shaded and unshaded leaves currentPatch%psn_z(CL,ft,iv) = 0._r8 ! psn is accumulated across sun and shaded leaves. @@ -624,9 +637,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) if (nint(c3psn(FT)) == 1)then - ci(cl,ft,iv) = 0.7_r8 * bc_in(s)%cair_pa(ifp) + ci(cl,ft,iv) = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) else - ci(cl,ft,iv) = 0.4_r8 * bc_in(s)%cair_pa(ifp) + ci(cl,ft,iv) = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) end if niter = 0 @@ -727,9 +740,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Convert gs_mol (umol H2O/m**2/s) to gs (m/s) and then to rs (s/m) gs = gs_mol / cf - if ( DEBUG ) write(iulog,*) 'EDPhoto 737 ', currentPatch%psn_z(cl,ft,iv) - if ( DEBUG ) write(iulog,*) 'EDPhoto 738 ', ag(cl,ft,iv) - if ( DEBUG ) write(iulog,*) 'EDPhoto 739 ', currentPatch%f_sun(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 737 ', currentPatch%psn_z(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', ag(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 739 ', currentPatch%f_sun(cl,ft,iv) !accumulate total photosynthesis umol/m2 ground/s-1. weight per unit sun and sha leaves. if(sunsha == 1)then !sunlit @@ -752,15 +765,15 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if - if ( DEBUG ) write(iulog,*) 'EDPhoto 758 ', currentPatch%psn_z(cl,ft,iv) - if ( DEBUG ) write(iulog,*) 'EDPhoto 759 ', ag(cl,ft,iv) - if ( DEBUG ) write(iulog,*) 'EDPhoto 760 ', currentPatch%f_sun(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', currentPatch%psn_z(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', ag(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 760 ', currentPatch%f_sun(cl,ft,iv) ! Make sure iterative solution is correct if (gs_mol < 0._r8) then - write (iulog,*)'Negative stomatal conductance:' - write (iulog,*)'ifp,iv,gs_mol= ',ifp,iv,gs_mol - call endrun(msg=errmsg(sourcefile, __LINE__)) + write (fates_log(),*)'Negative stomatal conductance:' + write (fates_log(),*)'ifp,iv,gs_mol= ',ifp,iv,gs_mol + call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Compare with Ball-Berry model: gs_mol = m * an * hs/cs p + b @@ -768,8 +781,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) gs_mol_err = bb_slope(ft)*max(an(cl,ft,iv), 0._r8)*hs/cs*bc_in(s)%forc_pbot + bbb(FT) if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then - write (iulog,*) 'CF: Ball-Berry error check - stomatal conductance error:' - write (iulog,*) gs_mol, gs_mol_err + write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' + write (fates_log(),*) gs_mol, gs_mol_err end if enddo !sunsha loop @@ -817,12 +830,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) !------------------------------------------------------------------------------ ! Convert from umolC/m2leaf/s to umolC/indiv/s ( x canopy area x 1m2 leaf area). tree_area = currentCohort%c_area/currentCohort%n - if ( DEBUG ) write(iulog,*) 'EDPhoto 816 ', currentCohort%gpp_tstep - if ( DEBUG ) write(iulog,*) 'EDPhoto 817 ', currentPatch%psn_z(cl,ft,1:currentCohort%nv-1) - if ( DEBUG ) write(iulog,*) 'EDPhoto 818 ', cl - if ( DEBUG ) write(iulog,*) 'EDPhoto 819 ', ft - if ( DEBUG ) write(iulog,*) 'EDPhoto 820 ', currentCohort%nv - if ( DEBUG ) write(iulog,*) 'EDPhoto 821 ', currentPatch%elai_profile(cl,ft,1:currentCohort%nv-1) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 816 ', currentCohort%gpp_tstep + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 817 ', currentPatch%psn_z(cl,ft,1:currentCohort%nv-1) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 818 ', cl + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 819 ', ft + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 820 ', currentCohort%nv + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 821 ', currentPatch%elai_profile(cl,ft,1:currentCohort%nv-1) if (currentCohort%nv > 1) then !is there canopy, and are the leaves on? @@ -844,7 +857,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if - if ( DEBUG ) write(iulog,*) 'EDPhoto 832 ', currentCohort%gpp_tstep + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 832 ', currentCohort%gpp_tstep laifrac = (currentCohort%treelai+currentCohort%treesai)-(currentCohort%nv-1)*dinc_ed @@ -852,18 +865,18 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%gscan = currentCohort%gscan+gs_cohort if ( DEBUG ) then - write(iulog,*) 'EDPhoto 868 ', currentCohort%gpp_tstep - write(iulog,*) 'EDPhoto 869 ', currentPatch%psn_z(cl,ft,currentCohort%nv) - write(iulog,*) 'EDPhoto 870 ', currentPatch%elai_profile(cl,ft,currentCohort%nv) - write(iulog,*) 'EDPhoto 871 ', laifrac - write(iulog,*) 'EDPhoto 872 ', tree_area - write(iulog,*) 'EDPhoto 873 ', currentCohort%nv, cl, ft + write(fates_log(),*) 'EDPhoto 868 ', currentCohort%gpp_tstep + write(fates_log(),*) 'EDPhoto 869 ', currentPatch%psn_z(cl,ft,currentCohort%nv) + write(fates_log(),*) 'EDPhoto 870 ', currentPatch%elai_profile(cl,ft,currentCohort%nv) + write(fates_log(),*) 'EDPhoto 871 ', laifrac + write(fates_log(),*) 'EDPhoto 872 ', tree_area + write(fates_log(),*) 'EDPhoto 873 ', currentCohort%nv, cl, ft endif currentCohort%gpp_tstep = currentCohort%gpp_tstep + currentPatch%psn_z(cl,ft,currentCohort%nv) * & currentPatch%elai_profile(cl,ft,currentCohort%nv) * laifrac * tree_area - if ( DEBUG ) write(iulog,*) 'EDPhoto 843 ', currentCohort%rdark + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 843 ', currentCohort%rdark currentCohort%rdark = currentCohort%rdark + lmr_z(cl,ft,currentCohort%nv) * & @@ -884,12 +897,15 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Calculate the amount of nitrogen in the above and below ground ! stem and root pools, used for maint resp + ! We are using the fine-root C:N ratio as an approximation for + ! the sapwood pools. + ! Units are in (kgN/plant) ! ------------------------------------------------------------------ live_agstem_n = ED_val_ag_biomass * currentCohort%bsw / & - pftcon%leafcn(currentCohort%pft) + frootcn(currentCohort%pft) live_bgstem_n = (1.0_r8-ED_val_ag_biomass) * currentCohort%bsw / & - pftcon%leafcn(currentCohort%pft) - froot_n = currentCohort%br / leafcn(currentCohort%pft) + frootcn(currentCohort%pft) + froot_n = currentCohort%br / frootcn(currentCohort%pft) !------------------------------------------------------------------------------ @@ -897,54 +913,48 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Leaf respn needs to be in the sub-layer loop to account for changing N through canopy. !------------------------------------------------------------------------------ - ! Above ground Live stem MR + ! Above ground Live stem MR (kgC/plant/s) ! ------------------------------------------------------------------ if (woody(ft) == 1) then tc = q10**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8) + ! kgC/s = kgN * kgC/kgN/s currentCohort%livestem_mr = live_agstem_n * base_mr_20 * tc - !convert from gC /indiv/s-1 to kgC/indiv/s-1 - currentCohort%livestem_mr = currentCohort%livestem_mr /1000.0_r8 else currentCohort%livestem_mr = 0._r8 end if - ! Fine Root MR + ! Fine Root MR (kgC/plant/s) ! ------------------------------------------------------------------ currentCohort%froot_mr = 0._r8 - do j = 1,nlevsoi + do j = 1,cp_numlevsoil tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%froot_mr = currentCohort%froot_mr + & froot_n * base_mr_20 * tcsoi * currentPatch%rootfr_ft(ft,j) enddo - ! convert from gC/indiv/s-1 to kgC/indiv/s-1 - currentCohort%froot_mr = currentCohort%froot_mr /1000.0_r8 - ! Coarse Root MR ! ------------------------------------------------------------------ if (woody(ft) == 1) then currentCohort%livecroot_mr = 0._r8 - do j = 1,nlevsoi + do j = 1,cp_numlevsoil ! Soil temperature used to adjust base rate of MR tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%livecroot_mr = currentCohort%livecroot_mr + & live_bgstem_n * base_mr_20 * tcsoi * & currentPatch%rootfr_ft(ft,j) enddo - ! convert from gC/indiv/s-1 to kgC/indiv/s-1 - currentCohort%livecroot_mr = currentCohort%livecroot_mr /1000.0_r8 else currentCohort%livecroot_mr = 0._r8 end if ! convert gpp from umol/indiv/s-1 to kgC/indiv/s-1 = X * 12 *10-6 * 10-3 - if ( DEBUG ) write(iulog,*) 'EDPhoto 904 ', currentCohort%resp_m - if ( DEBUG ) write(iulog,*) 'EDPhoto 905 ', currentCohort%rdark - if ( DEBUG ) write(iulog,*) 'EDPhoto 906 ', currentCohort%livestem_mr - if ( DEBUG ) write(iulog,*) 'EDPhoto 907 ', currentCohort%livecroot_mr - if ( DEBUG ) write(iulog,*) 'EDPhoto 908 ', currentCohort%froot_mr + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 904 ', currentCohort%resp_m + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 905 ', currentCohort%rdark + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 906 ', currentCohort%livestem_mr + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 907 ', currentCohort%livecroot_mr + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 908 ', currentCohort%froot_mr currentCohort%gpp_tstep = currentCohort%gpp_tstep * umolC_to_kgC ! add on whole plant respiration values in kgC/indiv/s-1 @@ -956,9 +966,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%resp_m = currentCohort%resp_m * dtime currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime - if ( DEBUG ) write(iulog,*) 'EDPhoto 911 ', currentCohort%gpp_tstep - if ( DEBUG ) write(iulog,*) 'EDPhoto 912 ', currentCohort%resp_tstep - if ( DEBUG ) write(iulog,*) 'EDPhoto 913 ', currentCohort%resp_m + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 911 ', currentCohort%gpp_tstep + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 912 ', currentCohort%resp_tstep + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 913 ', currentCohort%resp_m currentCohort%resp_g = ED_val_grperc(ft) * (max(0._r8,currentCohort%gpp_tstep - currentCohort%resp_m)) currentCohort%resp_tstep = currentCohort%resp_m + currentCohort%resp_g ! kgC/indiv/ts @@ -979,7 +989,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%gscan = 0._r8 end if else !pft<0 n<0 - write(iulog,*) 'CF: pft 0 or n 0',currentCohort%pft,currentCohort%n,currentCohort%indexnumber + write(fates_log(),*) 'CF: pft 0 or n 0',currentCohort%pft,currentCohort%n,currentCohort%indexnumber currentCohort%gpp_tstep = 0._r8 currentCohort%resp_m = 0._r8 currentCohort%gscan = 0._r8 @@ -1009,7 +1019,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if bc_out(s)%rssun_pa(ifp) = rscanopy bc_out(s)%rssha_pa(ifp) = rscanopy - bc_out(s)%gccanopy_pa(ifp) = 1.0_r8/rscanopy*cf/1000 !convert into umol m02 s-1 then mmol m-2 s-1. + bc_out(s)%gccanopy_pa(ifp) = 1.0_r8/rscanopy*cf/1000.0 !convert into umol m-2 s-1 then mmol m-2 s-1. end if currentPatch => currentPatch%younger @@ -1034,7 +1044,8 @@ function ft1_f(tl, ha) result(ans) ! 7/23/16: Copied over from CLM by Ryan Knox ! !!USES - use clm_varcon , only : rgas, tfrz + use FatesConstantsMod, only : rgas, & ! universal gas constant + tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm ! ! !ARGUMENTS: real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) @@ -1060,7 +1071,8 @@ function fth_f(tl,hd,se,scaleFactor) result(ans) ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 ! 7/23/16: Copied over from CLM by Ryan Knox ! - use clm_varcon , only : rgas, tfrz + use FatesConstantsMod, only : rgas, & ! universal gas constant + tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm ! ! !ARGUMENTS: real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) @@ -1088,8 +1100,10 @@ function fth25_f(hd,se)result(ans) ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 ! 7/23/16: Copied over from CLM by Ryan Knox ! - !!USES - use clm_varcon , only : rgas, tfrz + !!USES + use FatesConstantsMod, only : rgas, & ! universal gas constant + tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm + ! ! !ARGUMENTS: real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temperature function (J/mol) @@ -1131,8 +1145,8 @@ subroutine quadratic_f (a, b, c, r1, r2) !------------------------------------------------------------------------------ if (a == 0._r8) then - write (iulog,*) 'Quadratic solution error: a = ',a - call endrun(msg=errmsg(sourcefile, __LINE__)) + write (fates_log(),*) 'Quadratic solution error: a = ',a + call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (b >= 0._r8) then diff --git a/components/clm/src/ED/main/FatesConstantsMod.F90 b/components/clm/src/ED/main/FatesConstantsMod.F90 index 244f6f6505..8a9d6938e3 100644 --- a/components/clm/src/ED/main/FatesConstantsMod.F90 +++ b/components/clm/src/ED/main/FatesConstantsMod.F90 @@ -14,4 +14,33 @@ module FatesConstantsMod integer, parameter :: fates_short_string_length = 32 integer, parameter :: fates_long_string_length = 199 + + ! Unit conversion constants: + + ! Conversion factor umols of Carbon -> kg of Carbon (1 mol = 12g) + ! We do not use umolC_per_kg because it is a non-terminating decimal + real(fates_r8), parameter :: umolC_to_kgC = 12.0E-9_fates_r8 + + ! Conversion factor: grams per kilograms + real(fates_r8), parameter :: g_per_kg = 1000.0_fates_r8 + + ! Conversion factor: miligrams per grams + real(fates_r8), parameter :: mg_per_g = 1000.0_fates_r8 + + ! Conversion: secons per minute + real(fates_r8), parameter :: sec_per_min = 60.0_fates_r8 + + + ! Physical constants + + ! universal gas constant [J/K/kmole] + real(fates_r8), parameter :: rgas = 8314.4598_fates_r8 + + ! freezing point of water at 1 atm (K) + real(fates_r8), parameter :: t_water_freeze_k_1atm = 273.15_fates_r8 + + ! freezing point of water at triple point (K) + real(fates_r8), parameter :: t_water_freeze_k_triple = 273.16_fates_r8 + + end module FatesConstantsMod From fe9f7a11030653901411ce12e7f8f9235c69cb8e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 13 Oct 2016 18:29:11 -0700 Subject: [PATCH 08/11] intercellular co2 in photosynthesis was being represented as an array that was unnecessary, a scalar was sufficient. It was also calculated twice in the call sequence, where the first calculation was unnecessary and unused. --- .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 29 +++++++------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 79b4897728..6179e97fa5 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -102,8 +102,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: lmr_z(cp_nclmax,mxpft,cp_nlevcan) ! initial slope of CO2 response curve (C4 plants) real(r8) :: rs_z(cp_nclmax,mxpft,cp_nlevcan) ! stomatal resistance s/m real(r8) :: gs_z(cp_nclmax,mxpft,cp_nlevcan) ! stomatal conductance m/s - - real(r8) :: ci(cp_nclmax,mxpft,cp_nlevcan) ! intracellular leaf CO2 (Pa) + + real(r8) :: ci ! intracellular leaf CO2 (Pa) real(r8) :: lnc(mxpft) ! leaf N concentration (gN leaf/m^2) real(r8) :: kc( numpatchespercol ) ! Michaelis-Menten constant for CO2 (Pa) @@ -412,13 +412,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if bbb(FT) = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8) - ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 650 (RGK) - if (nint(c3psn(FT)) == 1)then - ci(:,FT,:) = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) - else - ci(:,FT,:) = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) - end if - ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) lnc(FT) = 1._r8 / (slatop(FT) * leafcn(FT)) @@ -637,9 +630,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) if (nint(c3psn(FT)) == 1)then - ci(cl,ft,iv) = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) + ci = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) else - ci(cl,ft,iv) = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) + ci = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) end if niter = 0 @@ -649,15 +642,15 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) niter = niter + 1 ! Save old ci - ciold = ci(cl,ft,iv) + ciold = ci ! Photosynthesis limitation rate calculations if (nint(c3psn(FT)) == 1)then ! C3: Rubisco-limited photosynthesis - ac = vcmax_z(cl,ft,iv) * max(ci(cl,ft,iv)-co2_cp(ifp), 0._r8) / (ci(cl,ft,iv)+kc(ifp)* & + ac = vcmax_z(cl,ft,iv) * max(ci-co2_cp(ifp), 0._r8) / (ci+kc(ifp)* & (1._r8+bc_in(s)%oair_pa(ifp)/ko(ifp))) ! C3: RuBP-limited photosynthesis - aj = je * max(ci(cl,ft,iv)-co2_cp(ifp), 0._r8) / (4._r8*ci(cl,ft,iv)+8._r8*co2_cp(ifp)) + aj = je * max(ci-co2_cp(ifp), 0._r8) / (4._r8*ci+8._r8*co2_cp(ifp)) ! C3: Product-limited photosynthesis ap = 3._r8 * tpu_z(cl,ft,iv) else @@ -681,7 +674,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if ! C4: PEP carboxylase-limited (CO2-limited) - ap = kp_z(cl,ft,iv) * max(ci(cl,ft,iv), 0._r8) / bc_in(s)%forc_pbot + ap = kp_z(cl,ft,iv) * max(ci, 0._r8) / bc_in(s)%forc_pbot end if ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap aquad = theta_cj(ps) @@ -715,14 +708,14 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) gs_mol = max(r1,r2) ! Derive new estimate for ci - ci(cl,ft,iv) = bc_in(s)%cair_pa(ifp) - an(cl,ft,iv) * bc_in(s)%forc_pbot * & + ci = bc_in(s)%cair_pa(ifp) - an(cl,ft,iv) * bc_in(s)%forc_pbot * & (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) ! Check for ci convergence. Delta ci/pair = mol/mol. Multiply by 10**6 to ! convert to umol/mol (ppm). Exit iteration if convergence criteria of +/- 1 x 10**-6 ppm ! is met OR if at least ten iterations (niter=10) are completed - if ((abs(ci(cl,ft,iv)-ciold)/bc_in(s)%forc_pbot*1.e06_r8 <= 2.e-06_r8) .or. niter == 5) then + if ((abs(ci-ciold)/bc_in(s)%forc_pbot*1.e06_r8 <= 2.e-06_r8) .or. niter == 5) then exitloop = 1 end if end do !iteration loop @@ -735,7 +728,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Final estimates for cs and ci (needed for early exit of ci iteration when an < 0) cs = bc_in(s)%cair_pa(ifp) - 1.4_r8/gb_mol * an(cl,ft,iv) * bc_in(s)%forc_pbot cs = max(cs,1.e-06_r8) - ci(cl,ft,iv) = bc_in(s)%cair_pa(ifp) - & + ci = bc_in(s)%cair_pa(ifp) - & an(cl,ft,iv) * bc_in(s)%forc_pbot * (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) ! Convert gs_mol (umol H2O/m**2/s) to gs (m/s) and then to rs (s/m) gs = gs_mol / cf From fb3e9ca8b79013d389af153a5c4ee5523bbafec3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 18 Oct 2016 13:33:47 -0700 Subject: [PATCH 09/11] renamed some local variables related to nitrogen content. --- .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 6179e97fa5..5e9c2ee39c 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -214,9 +214,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: rscanopy real(r8) :: elai - real(r8) :: live_agstem_n ! Live above-ground stem (sapwood) nitrogen content (gN/plant) - real(r8) :: live_bgstem_n ! Live below-ground stem (sapwood) nitrogen content (gN/plant) - real(r8) :: froot_n ! Fine root nitrogen content (gN/plant) + real(r8) :: live_stem_n ! Live stem (above-ground sapwood) nitrogen content (kgN/plant) + real(r8) :: live_croot_n ! Live coarse root (below-ground sapwood) nitrogen content (kgN/plant) + real(r8) :: froot_n ! Fine root nitrogen content (kgN/plant) ! Parameters ! ----------------------------------------------------------------------- @@ -894,9 +894,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! the sapwood pools. ! Units are in (kgN/plant) ! ------------------------------------------------------------------ - live_agstem_n = ED_val_ag_biomass * currentCohort%bsw / & + live_stem_n = ED_val_ag_biomass * currentCohort%bsw / & frootcn(currentCohort%pft) - live_bgstem_n = (1.0_r8-ED_val_ag_biomass) * currentCohort%bsw / & + live_croot_n = (1.0_r8-ED_val_ag_biomass) * currentCohort%bsw / & frootcn(currentCohort%pft) froot_n = currentCohort%br / frootcn(currentCohort%pft) @@ -906,12 +906,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Leaf respn needs to be in the sub-layer loop to account for changing N through canopy. !------------------------------------------------------------------------------ - ! Above ground Live stem MR (kgC/plant/s) + ! Live stem MR (kgC/plant/s) (above ground sapwood) ! ------------------------------------------------------------------ if (woody(ft) == 1) then tc = q10**((bc_in(s)%t_veg_pa(ifp)-tfrz - 20.0_r8)/10.0_r8) ! kgC/s = kgN * kgC/kgN/s - currentCohort%livestem_mr = live_agstem_n * base_mr_20 * tc + currentCohort%livestem_mr = live_stem_n * base_mr_20 * tc else currentCohort%livestem_mr = 0._r8 end if @@ -926,7 +926,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) froot_n * base_mr_20 * tcsoi * currentPatch%rootfr_ft(ft,j) enddo - ! Coarse Root MR + ! Coarse Root MR (kgC/plant/s) (below ground sapwood) ! ------------------------------------------------------------------ if (woody(ft) == 1) then currentCohort%livecroot_mr = 0._r8 @@ -934,7 +934,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Soil temperature used to adjust base rate of MR tcsoi = q10**((bc_in(s)%t_soisno_gl(j)-tfrz - 20.0_r8)/10.0_r8) currentCohort%livecroot_mr = currentCohort%livecroot_mr + & - live_bgstem_n * base_mr_20 * tcsoi * & + live_croot_n * base_mr_20 * tcsoi * & currentPatch%rootfr_ft(ft,j) enddo else From 44d762d950bb8c975af2bbe577cc857f076252ff Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 18 Oct 2016 14:47:02 -0700 Subject: [PATCH 10/11] Added cohort type variables size_class and size_by_pft_class. These variables were added to reduce the number of times these are re-calculated. They are in some cases needed for history diagnostics updated at the finest model time-scale on cohort variables, this was likely a very expensive calculation, so it was moved to be calculated only on the dynamics timestep. --- .../ED/biogeochem/EDCanopyStructureMod.F90 | 8 + .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 29 +++- components/clm/src/ED/main/EDTypesMod.F90 | 8 + components/clm/src/ED/main/HistoryIOMod.F90 | 138 +++++++++--------- 4 files changed, 111 insertions(+), 72 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 index 8ad6b71c79..c3ea17dd56 100755 --- a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 @@ -658,6 +658,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesInterfaceMod , only : bc_in_type use EDPatchDynamicsMod , only : set_patchno + use EDCohortDynamicsMod , only : size_and_type_class_index use EDGrowthFunctionsMod , only : tree_lai, c_area use EDEcophysConType , only : EDecophyscon use EDtypesMod , only : area @@ -709,6 +710,13 @@ subroutine canopy_summarization( nsites, sites, bc_in ) do while(associated(currentCohort)) ft = currentCohort%pft + + + ! Update the cohort's index within the size bin classes + ! Update the cohort's index within the SCPF classification system + call size_and_type_class_index(currentCohort%dbh,currentCohort%pft, & + currentCohort%size_class,currentCohort%size_by_pft_class) + currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore currentCohort%treelai = tree_lai(currentCohort) diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index d35afdb1d8..8f4ebd11f8 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -27,7 +27,7 @@ module EDCohortDynamicsMod public :: sort_cohorts public :: copy_cohort public :: count_cohorts -! public :: countCohorts + public :: size_and_type_class_index public :: allocate_live_biomass logical, parameter :: DEBUG = .false. ! local debug flag @@ -92,6 +92,9 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & new_cohort%balive = balive new_cohort%bstore = bstore + call size_and_type_class_index(new_cohort%dbh,new_cohort%pft, & + new_cohort%size_class,new_cohort%size_by_pft_class) + if ( DEBUG ) write(iulog,*) 'EDCohortDyn I ',bstore if (new_cohort%dbh <= 0.0_r8 .or. new_cohort%n == 0._r8 .or. new_cohort%pft == 0 & @@ -316,6 +319,8 @@ subroutine nan_cohort(cc_p) currentCohort%canopy_layer = 999 ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) currentCohort%NV = 999 ! Number of leaf layers: - currentCohort%status_coh = 999 ! growth status of plant (2 = leaves on , 1 = leaves off) + currentCohort%size_class = 999 ! size class index + currentCohort%size_by_pft_class = 999 ! size by pft classification index currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm @@ -1127,6 +1132,28 @@ function count_cohorts( currentPatch ) result ( backcount ) end function count_cohorts + ! ===================================================================================== + + subroutine size_and_type_class_index(dbh,pft,size_class,size_by_pft_class) + + use EDTypesMod, only: sclass_ed, & + nlevsclass_ed + + ! Arguments + real(r8),intent(in) :: dbh + integer,intent(in) :: pft + integer,intent(out) :: size_class + integer,intent(out) :: size_by_pft_class + + size_class = count(dbh-sclass_ed.ge.0.0) + + size_by_pft_class = (pft-1)*nlevsclass_ed+size_class + + return + end subroutine size_and_type_class_index + + + !-------------------------------------------------------------------------------------! ! function countCohorts( bounds, ed_allsites_inst ) result ( totNumCohorts ) ! diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 486db6150e..d02891cb28 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -189,6 +189,14 @@ module EDTypesMod real(r8) :: treesai ! stem area index of tree (total stem area (m2) / canopy area (m2) logical :: isnew ! flag to signify a new cohort, new cohorts have not experienced ! npp or mortality and should therefore not be fused or averaged + integer :: size_class ! An index that indicates which diameter size bin the cohort currently resides in + ! this is used for history output. We maintain this in the main cohort memory + ! because we don't want to continually re-calculate the cohort's position when + ! performing size diagnostics at high-frequency calls + integer :: size_by_pft_class ! An index that indicates the cohorts position of the joint size-class x functional + ! type classification. We also maintain this in the main cohort memory + ! because we don't want to continually re-calculate the cohort's position when + ! performing size diagnostics at high-frequency calls ! CARBON FLUXES real(r8) :: gpp ! GPP: kgC/indiv/year diff --git a/components/clm/src/ED/main/HistoryIOMod.F90 b/components/clm/src/ED/main/HistoryIOMod.F90 index e1a0b295f1..3b424d2edf 100644 --- a/components/clm/src/ED/main/HistoryIOMod.F90 +++ b/components/clm/src/ED/main/HistoryIOMod.F90 @@ -336,8 +336,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) integer :: lb1,ub1,lb2,ub2 ! IO array bounds for the calling thread integer :: ivar ! index of IO variable object vector integer :: ft ! functional type index - integer :: scpf ! index of the size-class x pft bin - integer :: sc ! index of the size-class bin real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column @@ -497,75 +495,76 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! ------------------------------------------------------------------------ dbh = ccohort%dbh !-0.5*(1./365.25)*ccohort%ddbhdt - sc = count(dbh-sclass_ed.ge.0.0) - scpf = (ft-1)*nlevsclass_ed+sc ! Flux Variables (cohorts must had experienced a day before any of these values ! have any meaning, otherwise they are just inialization values if( .not.(ccohort%isnew) ) then - hio_gpp_si_scpf(io_si,scpf) = hio_gpp_si_scpf(io_si,scpf) + & - n_perm2*ccohort%gpp ! [kgC/m2/yr] - hio_npp_totl_si_scpf(io_si,scpf) = hio_npp_totl_si_scpf(io_si,scpf) + & - ccohort%npp*n_perm2 - hio_npp_leaf_si_scpf(io_si,scpf) = hio_npp_leaf_si_scpf(io_si,scpf) + & - ccohort%npp_leaf*n_perm2 - hio_npp_fnrt_si_scpf(io_si,scpf) = hio_npp_fnrt_si_scpf(io_si,scpf) + & - ccohort%npp_froot*n_perm2 - hio_npp_bgsw_si_scpf(io_si,scpf) = hio_npp_bgsw_si_scpf(io_si,scpf) + & - ccohort%npp_bsw*(1._r8-ED_val_ag_biomass)*n_perm2 - hio_npp_agsw_si_scpf(io_si,scpf) = hio_npp_agsw_si_scpf(io_si,scpf) + & - ccohort%npp_bsw*ED_val_ag_biomass*n_perm2 - hio_npp_bgdw_si_scpf(io_si,scpf) = hio_npp_bgdw_si_scpf(io_si,scpf) + & - ccohort%npp_bdead*(1._r8-ED_val_ag_biomass)*n_perm2 - hio_npp_agdw_si_scpf(io_si,scpf) = hio_npp_agdw_si_scpf(io_si,scpf) + & - ccohort%npp_bdead*ED_val_ag_biomass*n_perm2 - hio_npp_seed_si_scpf(io_si,scpf) = hio_npp_seed_si_scpf(io_si,scpf) + & - ccohort%npp_bseed*n_perm2 - hio_npp_stor_si_scpf(io_si,scpf) = hio_npp_stor_si_scpf(io_si,scpf) + & - ccohort%npp_store*n_perm2 - - if( abs(ccohort%npp-(ccohort%npp_leaf+ccohort%npp_froot+ & - ccohort%npp_bsw+ccohort%npp_bdead+ & - ccohort%npp_bseed+ccohort%npp_store))>1.e-9) then - write(fates_log(),*) 'NPP Partitions are not balancing' - write(fates_log(),*) 'Fractional Error: ', & - abs(ccohort%npp-(ccohort%npp_leaf+ccohort%npp_froot+ & - ccohort%npp_bsw+ccohort%npp_bdead+ & - ccohort%npp_bseed+ccohort%npp_store))/ccohort%npp - write(fates_log(),*) 'Terms: ',ccohort%npp,ccohort%npp_leaf,ccohort%npp_froot, & - ccohort%npp_bsw,ccohort%npp_bdead, & - ccohort%npp_bseed,ccohort%npp_store - write(fates_log(),*) ' NPP components during FATES-HLM linking does not balance ' - stop ! we need termination control for FATES!!! - ! call endrun(msg=errMsg(__FILE__, __LINE__)) - end if + associate( scpf => ccohort%size_by_pft_class ) + + hio_gpp_si_scpf(io_si,scpf) = hio_gpp_si_scpf(io_si,scpf) + & + n_perm2*ccohort%gpp ! [kgC/m2/yr] + hio_npp_totl_si_scpf(io_si,scpf) = hio_npp_totl_si_scpf(io_si,scpf) + & + ccohort%npp*n_perm2 + hio_npp_leaf_si_scpf(io_si,scpf) = hio_npp_leaf_si_scpf(io_si,scpf) + & + ccohort%npp_leaf*n_perm2 + hio_npp_fnrt_si_scpf(io_si,scpf) = hio_npp_fnrt_si_scpf(io_si,scpf) + & + ccohort%npp_froot*n_perm2 + hio_npp_bgsw_si_scpf(io_si,scpf) = hio_npp_bgsw_si_scpf(io_si,scpf) + & + ccohort%npp_bsw*(1._r8-ED_val_ag_biomass)*n_perm2 + hio_npp_agsw_si_scpf(io_si,scpf) = hio_npp_agsw_si_scpf(io_si,scpf) + & + ccohort%npp_bsw*ED_val_ag_biomass*n_perm2 + hio_npp_bgdw_si_scpf(io_si,scpf) = hio_npp_bgdw_si_scpf(io_si,scpf) + & + ccohort%npp_bdead*(1._r8-ED_val_ag_biomass)*n_perm2 + hio_npp_agdw_si_scpf(io_si,scpf) = hio_npp_agdw_si_scpf(io_si,scpf) + & + ccohort%npp_bdead*ED_val_ag_biomass*n_perm2 + hio_npp_seed_si_scpf(io_si,scpf) = hio_npp_seed_si_scpf(io_si,scpf) + & + ccohort%npp_bseed*n_perm2 + hio_npp_stor_si_scpf(io_si,scpf) = hio_npp_stor_si_scpf(io_si,scpf) + & + ccohort%npp_store*n_perm2 + + if( abs(ccohort%npp-(ccohort%npp_leaf+ccohort%npp_froot+ & + ccohort%npp_bsw+ccohort%npp_bdead+ & + ccohort%npp_bseed+ccohort%npp_store))>1.e-9) then + write(fates_log(),*) 'NPP Partitions are not balancing' + write(fates_log(),*) 'Fractional Error: ', & + abs(ccohort%npp-(ccohort%npp_leaf+ccohort%npp_froot+ & + ccohort%npp_bsw+ccohort%npp_bdead+ & + ccohort%npp_bseed+ccohort%npp_store))/ccohort%npp + write(fates_log(),*) 'Terms: ',ccohort%npp,ccohort%npp_leaf,ccohort%npp_froot, & + ccohort%npp_bsw,ccohort%npp_bdead, & + ccohort%npp_bseed,ccohort%npp_store + write(fates_log(),*) ' NPP components during FATES-HLM linking does not balance ' + stop ! we need termination control for FATES!!! + ! call endrun(msg=errMsg(__FILE__, __LINE__)) + end if - ! Woody State Variables (basal area and number density and mortality) - if (pftcon%woody(ft) == 1) then - - hio_m1_si_scpf(io_si,scpf) = hio_m1_si_scpf(io_si,scpf) + ccohort%bmort*n_perm2*AREA - hio_m2_si_scpf(io_si,scpf) = hio_m2_si_scpf(io_si,scpf) + ccohort%hmort*n_perm2*AREA - hio_m3_si_scpf(io_si,scpf) = hio_m3_si_scpf(io_si,scpf) + ccohort%cmort*n_perm2*AREA - hio_m4_si_scpf(io_si,scpf) = hio_m4_si_scpf(io_si,scpf) + ccohort%imort*n_perm2*AREA - hio_m5_si_scpf(io_si,scpf) = hio_m5_si_scpf(io_si,scpf) + ccohort%fmort*n_perm2*AREA - - ! basal area [m2/ha] - hio_ba_si_scpf(io_si,scpf) = hio_ba_si_scpf(io_si,scpf) + & - 0.25_r8*3.14159_r8*((dbh/100.0_r8)**2.0_r8)*n_perm2*AREA - - ! number density [/ha] - hio_nplant_si_scpf(io_si,scpf) = hio_nplant_si_scpf(io_si,scpf) + AREA*n_perm2 - - ! Growth Incrments must have NaN check and woody check - if(ccohort%ddbhdt == ccohort%ddbhdt) then - hio_ddbh_si_scpf(io_si,scpf) = hio_ddbh_si_scpf(io_si,scpf) + & - ccohort%ddbhdt*n_perm2*AREA - else - hio_ddbh_si_scpf(io_si,scpf) = -999.9 - end if - end if - + ! Woody State Variables (basal area and number density and mortality) + if (pftcon%woody(ft) == 1) then + + hio_m1_si_scpf(io_si,scpf) = hio_m1_si_scpf(io_si,scpf) + ccohort%bmort*n_perm2*AREA + hio_m2_si_scpf(io_si,scpf) = hio_m2_si_scpf(io_si,scpf) + ccohort%hmort*n_perm2*AREA + hio_m3_si_scpf(io_si,scpf) = hio_m3_si_scpf(io_si,scpf) + ccohort%cmort*n_perm2*AREA + hio_m4_si_scpf(io_si,scpf) = hio_m4_si_scpf(io_si,scpf) + ccohort%imort*n_perm2*AREA + hio_m5_si_scpf(io_si,scpf) = hio_m5_si_scpf(io_si,scpf) + ccohort%fmort*n_perm2*AREA + + ! basal area [m2/ha] + hio_ba_si_scpf(io_si,scpf) = hio_ba_si_scpf(io_si,scpf) + & + 0.25_r8*3.14159_r8*((dbh/100.0_r8)**2.0_r8)*n_perm2*AREA + + ! number density [/ha] + hio_nplant_si_scpf(io_si,scpf) = hio_nplant_si_scpf(io_si,scpf) + AREA*n_perm2 + + ! Growth Incrments must have NaN check and woody check + if(ccohort%ddbhdt == ccohort%ddbhdt) then + hio_ddbh_si_scpf(io_si,scpf) = hio_ddbh_si_scpf(io_si,scpf) + & + ccohort%ddbhdt*n_perm2*AREA + else + hio_ddbh_si_scpf(io_si,scpf) = -999.9 + end if + end if + + end associate end if ccohort => ccohort%taller @@ -655,8 +654,6 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) integer :: lb1,ub1,lb2,ub2 ! IO array bounds for the calling thread integer :: ivar ! index of IO variable object vector integer :: ft ! functional type index - integer :: scpf ! index of the size-class x pft bin - integer :: sc ! index of the size-class bin real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column @@ -710,11 +707,9 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) endif if ( .not. ccohort%isnew ) then - + ! Calculate index for the scpf class - ft = ccohort%pft - sc = count(ccohort%dbh-sclass_ed.ge.0.0) - scpf = (ft-1)*nlevsclass_ed+sc + associate( scpf => ccohort%size_by_pft_class ) ! scale up cohort fluxes to their patches hio_npp_pa(io_pa) = hio_npp_pa(io_pa) + & @@ -761,6 +756,7 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + & ccohort%froot_mr * n_perm2 * daysecs * yeardays + end associate endif ccohort => ccohort%taller From 8e2b86f4e2383bced9476a2734d40365ec810cfc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 27 Oct 2016 15:56:34 -0700 Subject: [PATCH 11/11] Various syntactical fixes including: making rgas more explicit in its name, adding _r8 to local constants, adding an unset integer to global constants, removing some text that explained a constant. --- .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 22 +++--- .../src/ED/biogeophys/EDPhotosynthesisMod.F90 | 76 +++++++++---------- .../clm/src/ED/main/FatesConstantsMod.F90 | 18 ++++- 3 files changed, 66 insertions(+), 50 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 8f4ebd11f8..cdca9ec65b 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -293,6 +293,8 @@ subroutine nan_cohort(cc_p) ! ! !USES: use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) + use FatesConstantsMod, only : fates_unset_int + ! ! !ARGUMENTS type (ed_cohort_type), intent(inout), target :: cc_p @@ -314,13 +316,13 @@ subroutine nan_cohort(cc_p) nullify(currentCohort%siteptr) ! VEGETATION STRUCTURE - currentCohort%pft = 999 ! pft number - currentCohort%indexnumber = 999 ! unique number for each cohort. (within clump?) - currentCohort%canopy_layer = 999 ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) - currentCohort%NV = 999 ! Number of leaf layers: - - currentCohort%status_coh = 999 ! growth status of plant (2 = leaves on , 1 = leaves off) - currentCohort%size_class = 999 ! size class index - currentCohort%size_by_pft_class = 999 ! size by pft classification index + currentCohort%pft = fates_unset_int ! pft number + currentCohort%indexnumber = fates_unset_int ! unique number for each cohort. (within clump?) + currentCohort%canopy_layer = fates_unset_int ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) + currentCohort%NV = fates_unset_int ! Number of leaf layers: - + currentCohort%status_coh = fates_unset_int ! growth status of plant (2 = leaves on , 1 = leaves off) + currentCohort%size_class = fates_unset_int ! size class index + currentCohort%size_by_pft_class = fates_unset_int ! size by pft classification index currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm @@ -1136,8 +1138,8 @@ end function count_cohorts subroutine size_and_type_class_index(dbh,pft,size_class,size_by_pft_class) - use EDTypesMod, only: sclass_ed, & - nlevsclass_ed + use EDTypesMod, only: sclass_ed + use EDTypesMod, only: nlevsclass_ed ! Arguments real(r8),intent(in) :: dbh @@ -1145,7 +1147,7 @@ subroutine size_and_type_class_index(dbh,pft,size_class,size_by_pft_class) integer,intent(out) :: size_class integer,intent(out) :: size_by_pft_class - size_class = count(dbh-sclass_ed.ge.0.0) + size_class = count(dbh-sclass_ed.ge.0.0_r8) size_by_pft_class = (pft-1)*nlevsclass_ed+size_class diff --git a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 index 5e9c2ee39c..a9e6cf5049 100644 --- a/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 +++ b/components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90 @@ -45,50 +45,49 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! READS ARE REFACTORED (RGK 10-13-2016) use pftconMod , only : pftcon ! THIS WILL BE DEPRECATED WHEN PARAMETER ! READS ARE REFACTORED (RGK 10-13-2016) - use EDParamsMod , only : ED_val_grperc, & - ED_val_ag_biomass + use EDParamsMod , only : ED_val_grperc + use EDParamsMod , only : ED_val_ag_biomass use EDSharedParamsMod , only : EDParamsShareInst - use EDTypesMod , only : numpft_ed, & - dinc_ed, & - ed_patch_type, & - ed_cohort_type, & - ed_site_type, & - numpft_ed, & - numpatchespercol, & - cp_numlevsoil, & - cp_nlevcan, & - cp_nclmax + use EDTypesMod , only : numpft_ed + use EDTypesMod , only : dinc_ed + use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : ed_cohort_type + use EDTypesMod , only : ed_site_type + use EDTypesMod , only : numpft_ed + use EDTypesMod , only : numpatchespercol + use EDTypesMod , only : cp_numlevsoil + use EDTypesMod , only : cp_nlevcan + use EDTypesMod , only : cp_nclmax use EDEcophysContype , only : EDecophyscon - use FatesInterfaceMod , only : bc_in_type, & - bc_out_type + + use FatesInterfaceMod , only : bc_in_type + use FatesInterfaceMod , only : bc_out_type - use EDCanopyStructureMod,only: calc_areaindex + use EDCanopyStructureMod, only : calc_areaindex - use FatesConstantsMod, only : umolC_to_kgC, & ! micromole conversion to kgC - g_per_kg, & ! number of grams per kg - mg_per_g, & ! number of miligrams per g - sec_per_min, & ! seconds per minute (60!) - rgas, & ! universal gas constant - - tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atmosphere + use FatesConstantsMod, only : umolC_to_kgC + use FatesConstantsMod, only : g_per_kg + use FatesConstantsMod, only : mg_per_g + use FatesConstantsMod, only : sec_per_min + use FatesConstantsMod, only : umol_per_mmol + use FatesConstantsMod, only : rgas => rgas_J_K_kmol + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - ! ! !ARGUMENTS: + ! ----------------------------------------------------------------------------------- integer,intent(in) :: nsites type(ed_site_type),intent(inout),target :: sites(nsites) type(bc_in_type),intent(in) :: bc_in(nsites) type(bc_out_type),intent(inout) :: bc_out(nsites) real(r8),intent(in) :: dtime - ! - ! !CALLED FROM: - ! subroutine CanopyFluxes - ! + ! !LOCAL VARIABLES: + ! ----------------------------------------------------------------------------------- type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type), pointer :: currentCohort - ! + integer , parameter :: psn_type = 2 !c3 or c4. logical :: DEBUG = .false. @@ -237,8 +236,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! First guess on ratio between intracellular co2 and the atmosphere ! an iterator converges on actual - real(r8),parameter :: init_a2l_co2_c3 = 0.7 - real(r8),parameter :: init_a2l_co2_c4 = 0.4 + real(r8),parameter :: init_a2l_co2_c3 = 0.7_r8 + real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8 associate( & @@ -1012,7 +1011,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if bc_out(s)%rssun_pa(ifp) = rscanopy bc_out(s)%rssha_pa(ifp) = rscanopy - bc_out(s)%gccanopy_pa(ifp) = 1.0_r8/rscanopy*cf/1000.0 !convert into umol m-2 s-1 then mmol m-2 s-1. + bc_out(s)%gccanopy_pa(ifp) = 1.0_r8/rscanopy*cf/umol_per_mmol !convert into umol m-2 s-1 then mmol m-2 s-1. end if currentPatch => currentPatch%younger @@ -1037,8 +1036,8 @@ function ft1_f(tl, ha) result(ans) ! 7/23/16: Copied over from CLM by Ryan Knox ! !!USES - use FatesConstantsMod, only : rgas, & ! universal gas constant - tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm + use FatesConstantsMod, only : rgas => rgas_J_K_kmol + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm ! ! !ARGUMENTS: real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) @@ -1064,8 +1063,9 @@ function fth_f(tl,hd,se,scaleFactor) result(ans) ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 ! 7/23/16: Copied over from CLM by Ryan Knox ! - use FatesConstantsMod, only : rgas, & ! universal gas constant - tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm + use FatesConstantsMod, only : rgas => rgas_J_K_kmol + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + ! ! !ARGUMENTS: real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) @@ -1094,8 +1094,9 @@ function fth25_f(hd,se)result(ans) ! 7/23/16: Copied over from CLM by Ryan Knox ! !!USES - use FatesConstantsMod, only : rgas, & ! universal gas constant - tfrz => t_water_freeze_k_1atm ! Freezing point of water at 1 atm + + use FatesConstantsMod, only : rgas => rgas_J_K_kmol + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm ! ! !ARGUMENTS: @@ -1127,7 +1128,6 @@ subroutine quadratic_f (a, b, c, r1, r2) ! 7/23/16: Copied over from CLM by Ryan Knox ! ! !USES: - implicit none ! ! !ARGUMENTS: real(r8), intent(in) :: a,b,c ! Terms for quadratic equation diff --git a/components/clm/src/ED/main/FatesConstantsMod.F90 b/components/clm/src/ED/main/FatesConstantsMod.F90 index 8a9d6938e3..824d1ab2c3 100644 --- a/components/clm/src/ED/main/FatesConstantsMod.F90 +++ b/components/clm/src/ED/main/FatesConstantsMod.F90 @@ -15,6 +15,10 @@ module FatesConstantsMod integer, parameter :: fates_long_string_length = 199 + ! Unset and various other 'special' values + integer, parameter :: fates_unset_int = -9999 + + ! Unit conversion constants: ! Conversion factor umols of Carbon -> kg of Carbon (1 mol = 12g) @@ -27,14 +31,24 @@ module FatesConstantsMod ! Conversion factor: miligrams per grams real(fates_r8), parameter :: mg_per_g = 1000.0_fates_r8 + ! Conversion factor: micromoles per milimole + real(fates_r8), parameter :: umol_per_mmol = 1000.0_fates_r8 + + ! Conversion factor: milimoles per mole + real(fates_r8), parameter :: mmol_per_mol = 1000.0_fates_r8 + + ! Conversion factor: micromoles per mole + real(fates_r8), parameter :: umol_per_mol = 1.0E6_fates_r8 + + ! Conversion: secons per minute real(fates_r8), parameter :: sec_per_min = 60.0_fates_r8 ! Physical constants - ! universal gas constant [J/K/kmole] - real(fates_r8), parameter :: rgas = 8314.4598_fates_r8 + ! universal gas constant [J/K/kmol] + real(fates_r8), parameter :: rgas_J_K_kmol = 8314.4598_fates_r8 ! freezing point of water at 1 atm (K) real(fates_r8), parameter :: t_water_freeze_k_1atm = 273.15_fates_r8