From 6a78b7dab90a59ef82b133063b1e8e61eac933c1 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 1 Dec 2016 13:14:55 -0800 Subject: [PATCH 01/15] Partial refactors on photosynthesis. --- biogeophys/EDPhotosynthesisMod.F90 | 1018 +++++++++++++++++----------- 1 file changed, 624 insertions(+), 394 deletions(-) diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index a9e6cf5049..a0597d1ee5 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -105,9 +105,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) 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) - real(r8) :: ko( numpatchespercol ) ! Michaelis-Menten constant for O2 (Pa) - real(r8) :: co2_cp( numpatchespercol ) ! CO2 compensation point (Pa) + real(r8) :: kc ! Michaelis-Menten constant for CO2 (Pa) + real(r8) :: ko ! Michaelis-Menten constant for O2 (Pa) + real(r8) :: co2_cp ! CO2 compensation point (Pa) ! --------------------------------------------------------------- ! TO-DO: bbbopt is slated to be transferred to the parameter file @@ -127,17 +127,11 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) real(r8) :: kp25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C - real(r8) :: kc25 ! Michaelis-Menten constant for CO2 at 25C (Pa) - real(r8) :: ko25 ! Michaelis-Menten constant for O2 at 25C (Pa) - real(r8) :: cp25 ! CO2 compensation point at 25C (Pa) real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) real(r8) :: jmaxha ! activation energy for jmax (J/mol) real(r8) :: tpuha ! activation energy for tpu (J/mol) real(r8) :: lmrha ! activation energy for lmr (J/mol) - real(r8) :: kcha ! activation energy for kc (J/mol) - real(r8) :: koha ! activation energy for ko (J/mol) - real(r8) :: cpha ! activation energy for cp (J/mol) real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) @@ -172,7 +166,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: gs_mol ! leaf stomatal conductance (umol H2O/m**2/s) real(r8) :: gs ! leaf stomatal conductance (m/s) real(r8) :: hs ! fractional humidity at leaf surface (dimensionless) - real(r8) :: sco ! relative specificity of rubisco real(r8) :: tl ! leaf temperature in photosynthesis temperature function (K) real(r8) :: ha ! activation energy in photosynthesis temperature function (J/mol) real(r8) :: hd ! deactivation energy in photosynthesis temperature function (J/mol) @@ -193,9 +186,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) real(r8) :: ap ! product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m**2/s) - real(r8) :: ag(cp_nclmax,mxpft,cp_nlevcan) ! co-limited gross leaf photosynthesis (umol CO2/m**2/s) - 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. @@ -273,9 +264,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 ! except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282 - kcha = 79430._r8 - koha = 36380._r8 - cpha = 37830._r8 vcmaxha = 65330._r8 jmaxha = 43540._r8 tpuha = 53100._r8 @@ -316,101 +304,68 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) do s = 1,nsites + ! Multi-layer parameters scaled by leaf nitrogen profile. + ! Loop through each canopy layer to calculate nitrogen profile using + ! cumulative lai at the midpoint of the layer + ifp = 0 currentpatch => sites(s)%oldest_patch do while (associated(currentpatch)) - ifp = ifp+1 + ifp = ifp+1 + NCL_p = currentPatch%NCL_p + + ! Part I. Zero output boundary conditions + ! --------------------------------------------------------------------------- bc_out(s)%psncanopy_pa(ifp) = 0._r8 bc_out(s)%lmrcanopy_pa(ifp) = 0._r8 bc_out(s)%rssun_pa(ifp) = 0._r8 bc_out(s)%rssha_pa(ifp) = 0._r8 bc_out(s)%gccanopy_pa(ifp) = 0._r8 + ! Part II. Filter out patches ! Patch level filter flag for photosynthesis calculations ! has a short memory, flags: ! 1 = patch has not been called ! 2 = patch is currently marked for photosynthesis - ! 3 = patch has been called for photosynthesis at least once + ! 3 = patch has been called for photosynthesis already + ! --------------------------------------------------------------------------- if(bc_in(s)%filter_photo_pa(ifp)==2)then - currentPatch%ncan(:,:) = 0 - !redo the canopy structure algorithm to get round a bug that is happening for site 125, FT13. - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - - currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) - - currentCohort => currentCohort%shorter - - enddo !cohort - - currentPatch%nrad = currentPatch%ncan - do CL = 1,cp_nclmax - do ft = 1,numpft_ed - currentPatch%present(CL,ft) = 0 - do iv = 1, currentPatch%nrad(CL,ft); - if(currentPatch%canopy_area_profile(CL,ft,iv) > 0._r8)then - currentPatch%present(CL,ft) = 1 - end if - end do !iv - enddo !ft - enddo !CL - - ! kc, ko, currentPatch, from: Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 - ! - ! kc25 = 404.9 umol/mol - ! ko25 = 278.4 mmol/mol - ! cp25 = 42.75 umol/mol - ! - ! Derive sco from currentPatch and O2 using present-day O2 (0.209 mol/mol) and re-calculate - ! currentPatch to account for variation in O2 using currentPatch = 0.5 O2 / sco - ! + ! Part III. Calculate the number of sublayers for each pft and layer. And then identify + ! which layer/pft combinations have things in them. Output: + ! currentPatch%ncan(:,:) + ! currentPatch%present(:,:) + call UpdateCanopyNCanNRadPresent(currentPatch) + - kc25 = (404.9_r8 / 1.e06_r8) * bc_in(s)%forc_pbot - ko25 = (278.4_r8 / 1.e03_r8) * bc_in(s)%forc_pbot - sco = 0.5_r8 * 0.209_r8 / (42.75_r8 / 1.e06_r8) - cp25 = 0.5_r8 * bc_in(s)%oair_pa(ifp) / sco + ! Part IV. Identify some environmentally derived parameters: + ! Michaelis-Menten constant for CO2 (Pa) + ! Michaelis-Menten constant for O2 (Pa) + ! CO2 compensation point (Pa) + + call GetCanopyGasParameters(bc_in(s)%forc_pbot,bc_in(s)%oair_pa(ifp), & + bc_in(s)%t_veg_pa(ifp),kc,ko,co2_cp) + + + ! 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 - if(bc_in(s)%t_veg_pa(ifp).gt.150_r8.and.bc_in(s)%t_veg_pa(ifp).lt.350_r8)then - kc(ifp) = kc25 * ft1_f(bc_in(s)%t_veg_pa(ifp), kcha) - ko(ifp) = ko25 * ft1_f(bc_in(s)%t_veg_pa(ifp), koha) - co2_cp(ifp) = cp25 * ft1_f(bc_in(s)%t_veg_pa(ifp), cpha) - else - kc(ifp) = 1 - ko(ifp) = 1 - co2_cp(ifp) = 1 - end if + ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures + ! that hs does not go to zero. Also eair <= esat_tv so that hs <= 1 + ceair = min( max(bc_in(s)%eair_pa(ifp), 0.05_r8*bc_in(s)%esat_tv_pa(ifp)), bc_in(s)%esat_tv_pa(ifp) ) - end if - - currentpatch => currentpatch%younger - end do - ! Multi-layer parameters scaled by leaf nitrogen profile. - ! Loop through each canopy layer to calculate nitrogen profile using - ! cumulative lai at the midpoint of the layer - - ifp = 0 - currentpatch => sites(s)%oldest_patch - do while (associated(currentpatch)) - ifp = ifp+1 - if(bc_in(s)%filter_photo_pa(ifp)==2)then - - NCL_p = currentPatch%NCL_p - - do FT = 1,numpft_ed !calculate patch and pft specific propserties at canopy top. - - if (nint(c3psn(FT)) == 1)then - ps = 1 - else - ps = 2 - end if - bbb(FT) = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8) + ! Part V. Pre-process some variables that are PFT dependent + ! but not environmentally dependent + ! ------------------------------------------------------------------------ + do FT = 1,numpft_ed !calculate patch and pft specific properties at canopy top. + ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) lnc(FT) = 1._r8 / (slatop(FT) * leafcn(FT)) @@ -452,21 +407,20 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end do !FT + + ! If we are using plant hydro-dynamics, then several photosynthesis + ! variables will be available at the cohort scale, and not the + ! pft scale. So here we split and use different looping structures + ! ------------------------------------------------------------------ +! if ( use_fates_plant_hydro ) + + !==============================================================================! ! Calculate Nitrogen scaling factors and photosynthetic parameters. !==============================================================================! do CL = 1, NCL_p do FT = 1,numpft_ed - 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(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 - end if - enddo - if(currentPatch%present(CL,FT) == 1)then ! are there any leaves of this pft in this layer? if(CL==NCL_p)then !are we in the top canopy layer or a shaded layer? @@ -478,6 +432,15 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Loop through canopy layers (above snow). Respiration needs to be ! calculated every timestep. Others are calculated only if daytime do iv = 1, currentPatch%nrad(CL,FT) + + if (use_fates_plant_hydro) then + !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) + !! btran = currentCohort%btran(iv) + else + !! bbb = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8) + !! btran = currentPatch%btran_ft(currentCohort%pft) + end if + vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) !vegetation area index. if (iv == 1) then laican = laican + 0.5_r8 * vai @@ -485,313 +448,45 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) laican = laican + 0.5_r8 * (currentPatch%elai_profile(CL,FT,iv-1)+ & currentPatch%esai_profile(CL,FT,iv-1))+vai end if - + ! Scale for leaf nitrogen profile nscaler = exp(-kn(FT) * laican) - ! Maintenance respiration: umol CO2 / m**2 [leaf] / s - lmr25 = lmr25top(FT) * nscaler - - if (nint(c3psn(FT)) == 1)then - lmr_z(CL,FT,iv) = lmr25 * ft1_f(bc_in(s)%t_veg_pa(ifp), lmrha) * & - fth_f(bc_in(s)%t_veg_pa(ifp), lmrhd, lmrse, lmrc) - else - lmr_z(CL,FT,iv) = lmr25 * 2._r8**((bc_in(s)%t_veg_pa(ifp)-(tfrz+25._r8))/10._r8) - lmr_z(CL,FT,iv) = lmr_z(CL,FT,iv) / (1._r8 + exp( 1.3_r8*(bc_in(s)%t_veg_pa(ifp)-(tfrz+55._r8)) )) - end if - - if (currentPatch%ed_parsun_z(CL,FT,iv) <= 0._r8) then ! night time - vcmax_z(CL,FT,iv) = 0._r8 - jmax_z(CL,FT,iv) = 0._r8 - tpu_z(CL,FT,iv) = 0._r8 - kp_z(CL,FT,iv) = 0._r8 - else ! day time - vcmax25 = vcmax25top(FT) * nscaler - jmax25 = jmax25top(FT) * nscaler - tpu25 = tpu25top(FT) * nscaler - kp25 = kp25top(FT) * nscaler - - ! Adjust for temperature - vcmax_z(CL,FT,iv) = vcmax25 * ft1_f(bc_in(s)%t_veg_pa(ifp), vcmaxha) * & - fth_f(bc_in(s)%t_veg_pa(ifp), vcmaxhd, vcmaxse, vcmaxc) - jmax_z(CL,FT,iv) = jmax25 * ft1_f(bc_in(s)%t_veg_pa(ifp), jmaxha) * & - fth_f(bc_in(s)%t_veg_pa(ifp), jmaxhd, jmaxse, jmaxc) - tpu_z(CL,FT,iv) = tpu25 * ft1_f(bc_in(s)%t_veg_pa(ifp), tpuha) * & - fth_f(bc_in(s)%t_veg_pa(ifp), tpuhd, tpuse, tpuc) - - if (nint(c3psn(FT)) /= 1) then - vcmax_z(CL,FT,iv) = vcmax25 * 2._r8**((bc_in(s)%t_veg_pa(ifp)-(tfrz+25._r8))/10._r8) - vcmax_z(CL,FT,iv) = vcmax_z(CL,FT,iv) / (1._r8 + & - exp( 0.2_r8*((tfrz+15._r8)-bc_in(s)%t_veg_pa(ifp)) )) - vcmax_z(CL,FT,iv) = vcmax_z(CL,FT,iv) / (1._r8 + & - exp( 0.3_r8*(bc_in(s)%t_veg_pa(ifp)-(tfrz+40._r8)) )) - end if - kp_z(CL,FT,iv) = kp25 * 2._r8**((bc_in(s)%t_veg_pa(ifp)-(tfrz+25._r8))/10._r8) !q10 response of product limited psn. - end if - ! Adjust for soil water:(umol co2/m**2/s) + call LeafLayerPhotosynthesis(currentPatch%ed_parsun_z(CL,FT,iv), & ! in + currentPatch%ed_parsha_z(CL,FT,iv), & ! in + currentPatch%ed_laisun_z(CL,FT,iv), & ! in + currentPatch%ed_laisha_z(CL,FT,iv), & ! in + currentPatch%canopy_area_profile(CL,FT,iv), & ! in + ft, & ! in + nscaler, & ! in + lmr25top(ft), & ! in + vcmax25top(ft), & ! in + jmax25top(ft), & ! in + tpu25top(ft), & ! in + kp25top(ft), & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + btran, & ! in + bbb, & ! in + cf, & ! in + gb_mol, & ! in + lmr_z(CL,FT,iv), & ! out + currentPatch%psn_z(cl,ft,iv), & ! out + rs_z(CL,FT,iv), & ! out + an_av(CL,FT,iv)) ! out + - vcmax_z(CL,FT,iv) = vcmax_z(CL,FT,iv) * currentPatch%btran_ft(FT) - ! completely removed respiration drought response - ! - (lmr_z(CL,FT,iv) * (1.0_r8-currentPatch%btran_ft(FT)) *pftcon%resp_drought_response(FT)) - lmr_z(CL,FT,iv) = lmr_z(CL,FT,iv) + + end do ! iv end if !present enddo !PFT enddo !CL - !==============================================================================! - ! Leaf-level photosynthesis and stomatal conductance - !==============================================================================! - ! 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 - - ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures - ! that hs does not go to zero. Also eair <= esat_tv so that hs <= 1 - ceair = min( max(bc_in(s)%eair_pa(ifp), 0.05_r8*bc_in(s)%esat_tv_pa(ifp)), bc_in(s)%esat_tv_pa(ifp) ) - - ! Loop through canopy layers (above snow). Only do calculations if daytime - do CL = 1, NCL_p - do FT = 1,numpft_ed - if (nint(c3psn(FT)) == 1)then - ps = 1 - else - ps = 2 - 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(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 - aj = 0._r8 - ap = 0._r8 - ag(CL,FT,iv) = 0._r8 - an(CL,FT,iv) = ag(CL,FT,iv) - lmr_z(CL,FT,iv) - an_av(cl,ft,iv) = 0._r8 - currentPatch%psn_z(cl,ft,iv) = 0._r8 - rs_z(CL,FT,iv) = min(rsmax0, 1._r8/bbb(FT) * cf) - - else ! day time - !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) - - 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(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. - rs_z(CL,FT,iv) = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. - gs_z(CL,FT,iv) = 0._r8 - an_av(CL,FT,iv) = 0._r8 - do sunsha = 1,2 - ! Electron transport rate for C3 plants. Convert par from W/m2 to umol photons/m**2/s - ! using the factor 4.6 - ! Convert from units of par absorbed per unit ground area to par absorbed per unit leaf area. - - if(sunsha == 1)then !sunlit - if((currentPatch%ed_laisun_z(CL,FT,iv) * currentPatch%canopy_area_profile(CL,FT,iv)) > & - 0.0000000001_r8)then - - qabs = currentPatch%ed_parsun_z(CL,FT,iv) / (currentPatch%ed_laisun_z(CL,FT,iv) * & - currentPatch%canopy_area_profile(CL,FT,iv)) - qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 - - else - qabs = 0.0_r8 - end if - else - - qabs = currentPatch%ed_parsha_z(CL,FT,iv) / (currentPatch%ed_laisha_z(CL,FT,iv) * & - currentPatch%canopy_area_profile(CL,FT,iv)) - qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 - - end if - - !convert the absorbed par into absorbed par per m2 of leaf, - ! so it is consistant with the vcmax and lmr numbers. - aquad = theta_psii - bquad = -(qabs + jmax_z(cl,ft,iv)) - cquad = qabs * jmax_z(cl,ft,iv) - call quadratic_f (aquad, bquad, cquad, r1, r2) - je = min(r1,r2) - - ! Iterative loop for ci beginning with initial guess - ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) - - if (nint(c3psn(FT)) == 1)then - ci = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) - else - ci = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) - end if - - niter = 0 - exitloop = 0 - do while(exitloop == 0) - ! Increment iteration counter. Stop if too many iterations - niter = niter + 1 - - ! Save old ci - ciold = ci - - ! Photosynthesis limitation rate calculations - if (nint(c3psn(FT)) == 1)then - ! C3: Rubisco-limited photosynthesis - 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-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 - ! C4: Rubisco-limited photosynthesis - ac = vcmax_z(cl,ft,iv) - ! C4: RuBP-limited photosynthesis - if(sunsha == 1)then !sunlit - if((currentPatch%ed_laisun_z(cl,ft,iv) * currentPatch%canopy_area_profile(cl,ft,iv)) > & - 0.0000000001_r8)then !guard against /0's in the night. - aj = qe(ps) * currentPatch%ed_parsun_z(cl,ft,iv) * 4.6_r8 - !convert from per cohort to per m2 of leaf) - aj = aj / (currentPatch%ed_laisun_z(cl,ft,iv) * & - currentPatch%canopy_area_profile(cl,ft,iv)) - else - aj = 0._r8 - end if - else - aj = qe(ps) * currentPatch%ed_parsha_z(cl,ft,iv) * 4.6_r8 - aj = aj / (currentPatch%ed_laisha_z(cl,ft,iv) * & - currentPatch%canopy_area_profile(cl,ft,iv)) - end if - - ! C4: PEP carboxylase-limited (CO2-limited) - 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) - bquad = -(ac + aj) - cquad = ac * aj - call quadratic_f (aquad, bquad, cquad, r1, r2) - ai = min(r1,r2) - - aquad = theta_ip - bquad = -(ai + ap) - cquad = ai * ap - call quadratic_f (aquad, bquad, cquad, r1, r2) - ag(cl,ft,iv) = min(r1,r2) - - ! Net carbon assimilation. Exit iteration if an < 0 - an(cl,ft,iv) = ag(cl,ft,iv) - lmr_z(cl,ft,iv) - if (an(cl,ft,iv) < 0._r8) then - exitloop = 1 - end if - - ! Quadratic gs_mol calculation with an known. Valid for an >= 0. - ! With an <= 0, then gs_mol = bbb - - 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) - aquad = cs - bquad = cs*(gb_mol - bbb(FT)) - bb_slope(ft)*an(cl,ft,iv)*bc_in(s)%forc_pbot - cquad = -gb_mol*(cs*bbb(FT) + & - bb_slope(ft)*an(cl,ft,iv)*bc_in(s)%forc_pbot*ceair/bc_in(s)%esat_tv_pa(ifp)) - call quadratic_f (aquad, bquad, cquad, r1, r2) - gs_mol = max(r1,r2) - - ! Derive new estimate for ci - 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-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 - - ! End of ci iteration. Check for an < 0, in which case gs_mol = bbb - if (an(cl,ft,iv) < 0._r8) then - gs_mol = bbb(FT) - end if - - ! 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 = 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 - - 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 - - currentPatch%psn_z(cl,ft,iv) = currentPatch%psn_z(cl,ft,iv) + ag(cl,ft,iv) * & - currentPatch%f_sun(cl,ft,iv) - an_av(cl,ft,iv) = an_av(cl,ft,iv) + an(cl,ft,iv) * & - currentPatch%f_sun(cl,ft,iv) - gs_z(cl,ft,iv) = gs_z(cl,ft,iv) + 1._r8/(min(1._r8/gs, rsmax0)) * & - currentPatch%f_sun(cl,ft,iv) - - else - - currentPatch%psn_z(cl,ft,iv) = currentPatch%psn_z(cl,ft,iv) + ag(cl,ft,iv) & - * (1.0_r8-currentPatch%f_sun(cl,ft,iv)) - an_av(cl,ft,iv) = an_av(cl,ft,iv) + an(cl,ft,iv) & - * (1.0_r8-currentPatch%f_sun(cl,ft,iv)) - gs_z(cl,ft,iv) = gs_z(cl,ft,iv) + & - 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-currentPatch%f_sun(cl,ft,iv)) - - end if - - 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 (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 - hs = (gb_mol*ceair + gs_mol*bc_in(s)%esat_tv_pa(ifp)) / ((gb_mol+gs_mol)*bc_in(s)%esat_tv_pa(ifp)) - 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 (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' - write (fates_log(),*) gs_mol, gs_mol_err - end if - - enddo !sunsha loop - !average leaf-level stomatal resistance rate over sun and shade leaves... - rs_z(cl,ft,iv) = 1._r8/gs_z(cl,ft,iv) - else !No leaf area. This layer is present only because of stems. (leaves are off, or have reduced to 0 - currentPatch%psn_z(cl,ft,iv) = 0._r8 - rs_z(CL,FT,iv) = min(rsmax0, 1._r8/bbb(FT) * cf) - - end if !is there leaf area? - - - end if ! night or day - end do ! iv canopy layer - end if ! present(L,ft) ? rd_array - end do ! PFT loop - end do !canopy layer + !==============================================================================! ! Unpack fluxes from arrays into cohorts @@ -883,9 +578,17 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! THIS CALCULATION SHOULD BE MOVED TO THE ALLOMETRY MODULE (RGK 10-8-2016) + ! ------ IT ALSO SHOULD ALREADY HAVE BEEN CALCULATED RIGHT? + ! ------ CHANGING TO A CHECK ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - currentCohort%bsw = EDecophyscon%sapwood_ratio(currentCohort%pft) * & - currentCohort%hite * (currentCohort%balive + currentCohort%laimemory)*leaf_frac + if ( abs(currentCohort%bsw - (EDecophyscon%sapwood_ratio(currentCohort%pft) * & + currentCohort%hite * (currentCohort%balive + currentCohort%laimemory)*leaf_frac) ) & + > 1e-9 ) then + write(fates_log(),*) 'Sapwood biomass calculated during photosynthesis' + write(fates_log(),*) 'does not match what is contained in cohort%bsw' + write(fates_log(),*) 'which is the prognostic variable. Stopping.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Calculate the amount of nitrogen in the above and below ground ! stem and root pools, used for maint resp @@ -1021,11 +724,396 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end do !site loop end associate - end subroutine Photosynthesis_ED ! ======================================================================================= +subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in + parsha_lsl, & ! in + laisun_lsl, & ! in + laisha_lsl, & ! in + canopy_area_lsl, & ! in + ft, & ! in + nscaler, & ! in + lmr25top_ft, & ! in + vcmax25top_ft, & ! in + jmax25top_ft, & ! in + tpu25top_ft, & ! in + kp25top_ft, & ! in + t_veg, & ! in + btran, & ! in + bbb, & ! in + cf, & ! in + gb_mol, & ! in + lmr_out, & ! out + psn_out, & ! out + rstoma_out, & ! out + anet_av_out) ! out + + ! ------------------------------------------------------------------------------------ + ! This subroutine calculates photosynthesis and stomatal conductance within each leaf + ! sublayer. + ! A note on naming conventions: As this subroutine is called for every + ! leaf-sublayer, many of the arguments are specific to that "leaf sub layer" + ! (LSL), those variables are given a dimension tag "_lsl" + ! Other arguments or variables may be indicative of scales broader than the LSL. + ! ------------------------------------------------------------------------------------ + + ! Arguments + ! ------------------------------------------------------------------------ + real(r8), intent(in) :: parsun_lsl + real(r8), intent(in) :: parsha_lsl + real(r8), intent(in) :: laisun_lsl + real(r8), intent(in) :: laisha_lsl + real(r8), intent(in) :: elai_lsl + real(r8), intent(in) :: esai_lsl + real(r8), intent(in) :: canopy_area_lsl + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: lmr25top_ft ! canopy top leaf maint resp rate at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C + ! for this pft (umol electrons/m**2/s) + real(r8), intent(in) :: tpu25top_ft ! canopy top triose phosphate utilization rate at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve + ! (C4 plants) at 25C, canopy top, this pft + real(r8), intent(in) :: t_veg ! vegetation temperature + real(r8), intent(in) :: btran ! + real(r8), intent(in) :: bbb + real(r8), intent(in) :: cf + real(r8), intent(in) :: gb_mol + + real(r8), intent(out) :: lmr_out + real(r8), intent(out) :: psn_out + real(r8), intent(out) :: rstoma_out ! stomatal resistance (1/gs_lsl) (s/m) + real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) averaged over sun and shade leaves. + real(r8), intent(out) :: gstoma_out ! Stomatal Conductance of this leaf layer (m/s) + + + + ! Locals + ! ------------------------------------------------------------------------ + integer :: ps ! Index for the different photosynthetic pathways C3,C4 + integer :: sunsha ! Index for differentiating sun and shade + + real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C (umol CO2/m**2/s) + real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C (umol electrons/m**2/s) + real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) + real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) + real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C + + + !!!!!!! NOTE: WE CAN REMOVE THE LSL IDENTIFIER FROM LOCALS AFTER THEY HAVE MORE VERBOSE NAMES + ! + + real(r8) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) + real(r8) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) + + + real(r8) :: agross ! co-limited gross leaf photosynthesis (umol CO2/m**2/s) + real(r8) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) + + + + ! Parameters + ! ------------------------------------------------------------------------ + + associate( & + c3psn => pftcon%c3psn , & ! photosynthetic pathway: 0. = c4, 1. = c3 + slatop => pftcon%slatop , & ! specific leaf area at top of canopy, projected area basis [m^2/gC] + flnr => pftcon%flnr , & ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) + 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 + + + if (nint(c3psn(ft)) == 1)then + ps = 1 + else + ps = 2 + end if + + + + ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s + ! ---------------------------------------------------------------------------------- + lmr25 = lmr25top_ft * nscaler + if ( nint(c3psn(ft)) == 1)then + lmr_out = lmr25 * ft1_f(t_veg, lmrha) * & + fth_f(t_veg, lmrhd, lmrse, lmrc) + else + lmr_out = lmr25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) + lmr_out = lmr_out / (1._r8 + exp( 1.3_r8*(t_veg-(tfrz+55._r8)) )) + end if + + + + + ! Part II: Localized Biophysical Rates + ! ---------------------------------------------------------------------------------- + + if ( parsun_lsl <= 0._r8) then ! night time + vcmax = 0._r8 + jmax = 0._r8 + tpu = 0._r8 + co2_rcurve_islope = 0._r8 + else ! day time + vcmax25 = vcmax25top_ft * nscaler + jmax25 = jmax25top_ft * nscaler + tpu25 = tpu25top_ft * nscaler + co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler + + ! Adjust for temperature + vcmax = vcmax25 * ft1_f(t_veg, vcmaxha) * fth_f(t_veg, vcmaxhd, vcmaxse, vcmaxc) + jmax = jmax25 * ft1_f(t_veg, jmaxha) * fth_f(t_veg, jmaxhd, jmaxse, jmaxc) + tpu = tpu25 * ft1_f(t_veg, tpuha) * fth_f(t_veg, tpuhd, tpuse, tpuc) + + if (nint(c3psn(FT)) /= 1) then + vcmax = vcmax25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) + vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-t_veg ) )) + vcmax = vcmax / (1._r8 + exp( 0.3_r8*(t_veg-(tfrz+40._r8)) )) + end if + co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) !q10 response of product limited psn. + end if + + ! Adjust for water limitations + vcmax = vcmax * btran + + ! Leaf Maintenance Respiration has no direct water limitation effect + ! lmr_out = lmr_out * (nothing) + + + + + + + ! Part III: Photosynthesis and Conductance + ! ---------------------------------------------------------------------------------- + + if ( parsun_lsl <= 0._r8 ) then ! night time + + anet_av_out = 0._r8 + psn_out = 0._r8 + rstoma_out = min(rsmax0, 1._r8/bbb * cf) + + else ! day time (a little bit more complicated ...) + + if ( DEBUG ) write(fates_log(),*) 'EDphot 594 ',laisun_lsl + if ( DEBUG ) write(fates_log(),*) 'EDphot 595 ',laisha_lsl + + !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) + if ( laisun_lsl + laisha_lsl > 0._r8 ) then + + if ( DEBUG ) write(fates_log(),*) '600 in laisun, laisha loop ' + + !Loop aroun shaded and unshaded leaves + psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. + rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. + anet_av_out = 0._r8 + gstoma_out = 0._r8 + + do sunsha = 1,2 + ! Electron transport rate for C3 plants. Convert par from W/m2 to umol photons/m**2/s + ! using the factor 4.6 + ! Convert from units of par absorbed per unit ground area to par absorbed per unit leaf area. + + if(sunsha == 1)then !sunlit + if(( laisun_lsl * canopy_area_lsl) > 0.0000000001_r8)then + + qabs = parsun_lsl / (laisun_lsl * canopy_area_profile_lsl ) + qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 + + else + qabs = 0.0_r8 + end if + else + + qabs = currentPatch%ed_parsha_lsl / (currentPatch%ed_laisha_lsl * & + currentPatch%canopy_area_profile_lsl) + qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 + + end if + + !convert the absorbed par into absorbed par per m2 of leaf, + ! so it is consistant with the vcmax and lmr numbers. + aquad = theta_psii + bquad = -(qabs + jmax_z(cl,ft,iv)) + cquad = qabs * jmax_z(cl,ft,iv) + call quadratic_f (aquad, bquad, cquad, r1, r2) + je = min(r1,r2) + + ! Iterative loop for ci beginning with initial guess + ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) + + if (nint(c3psn(FT)) == 1)then + ci = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) + else + ci = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) + end if + + niter = 0 + exitloop = 0 + do while(exitloop == 0) + ! Increment iteration counter. Stop if too many iterations + niter = niter + 1 + + ! Save old ci + ciold = ci + + ! Photosynthesis limitation rate calculations + if (nint(c3psn(FT)) == 1)then + ! C3: Rubisco-limited photosynthesis + ac = vcmax_z(cl,ft,iv) * max(ci-co2_cp, 0._r8) / (ci+kc * & + (1._r8+bc_in(s)%oair_pa(ifp)/ko)) + ! C3: RuBP-limited photosynthesis + aj = je * max(ci-co2_cp, 0._r8) / (4._r8*ci+8._r8*co2_cp) + ! C3: Product-limited photosynthesis + ap = 3._r8 * tpu_z(cl,ft,iv) + else + ! C4: Rubisco-limited photosynthesis + ac = vcmax_z(cl,ft,iv) + ! C4: RuBP-limited photosynthesis + if(sunsha == 1)then !sunlit + if((currentPatch%ed_laisun_z(cl,ft,iv) * currentPatch%canopy_area_profile(cl,ft,iv)) > & + 0.0000000001_r8)then !guard against /0's in the night. + aj = qe(ps) * currentPatch%ed_parsun_z(cl,ft,iv) * 4.6_r8 + !convert from per cohort to per m2 of leaf) + aj = aj / (currentPatch%ed_laisun_z(cl,ft,iv) * & + currentPatch%canopy_area_profile(cl,ft,iv)) + else + aj = 0._r8 + end if + else + aj = qe(ps) * currentPatch%ed_parsha_z(cl,ft,iv) * 4.6_r8 + aj = aj / (currentPatch%ed_laisha_z(cl,ft,iv) * & + currentPatch%canopy_area_profile(cl,ft,iv)) + end if + + ! C4: PEP carboxylase-limited (CO2-limited) + ap = co2_rcurve_islope_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) + bquad = -(ac + aj) + cquad = ac * aj + call quadratic_f (aquad, bquad, cquad, r1, r2) + ai = min(r1,r2) + + aquad = theta_ip + bquad = -(ai + ap) + cquad = ai * ap + call quadratic_f (aquad, bquad, cquad, r1, r2) + agross = min(r1,r2) + + ! Net carbon assimilation. Exit iteration if an < 0 + anet = agross - lmr_out + if (anet < 0._r8) then + exitloop = 1 + end if + + ! Quadratic gs_mol calculation with an known. Valid for an >= 0. + ! With an <= 0, then gs_mol = bbb + + cs = bc_in(s)%cair_pa(ifp) - 1.4_r8/gb_mol * anet * bc_in(s)%forc_pbot + cs = max(cs,1.e-06_r8) + aquad = cs + bquad = cs*(gb_mol - bbb(FT)) - bb_slope(ft)*anet*bc_in(s)%forc_pbot + cquad = -gb_mol*(cs*bbb(FT) + & + bb_slope(ft)*anet*bc_in(s)%forc_pbot*ceair/bc_in(s)%esat_tv_pa(ifp)) + call quadratic_f (aquad, bquad, cquad, r1, r2) + gs_mol = max(r1,r2) + + ! Derive new estimate for ci + ci = bc_in(s)%cair_pa(ifp) - anet * 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-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 + + ! End of ci iteration. Check for an < 0, in which case gs_mol = bbb + if (anet < 0._r8) then + gs_mol = bbb(FT) + end if + + ! 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 * anet * bc_in(s)%forc_pbot + cs = max(cs,1.e-06_r8) + ci = bc_in(s)%cair_pa(ifp) - & + anet * 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 + + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 737 ', currentPatch%psn_z(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross + 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 + + psn_out = psn_out + agross * f_sun_lsl + anet_av_out = anet_av_out + anet * f_sun_lsl + gstoma_out = gstoma_out + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl + + else + + psn_out = psn_out + agross * (1.0_r8-f_sun_lsl) + anet_av_out = anet_av + anet * (1.0_r8-f_sun_lsl) + gstoma_out = gstoma_out + & + 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) + + end if + + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', currentPatch%psn_z(cl,ft,iv) + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', agross + 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 (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 + hs = (gb_mol*ceair + gs_mol*bc_in(s)%esat_tv_pa(ifp)) / ((gb_mol+gs_mol)*bc_in(s)%esat_tv_pa(ifp)) + gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/cs*bc_in(s)%forc_pbot + bbb(FT) + + if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then + write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' + write (fates_log(),*) gs_mol, gs_mol_err + end if + + enddo !sunsha loop + + !average leaf-level stomatal resistance rate over sun and shade leaves... + rstoma_out = 1._r8/gstoma_out + + else !No leaf area. This layer is present only because of stems. (leaves are off, or have reduced to 0 + currentPatch%psn_z(cl,ft,iv) = 0._r8 + rs_lsl = min(rsmax0, 1._r8/bbb(FT) * cf) + + end if !is there leaf area? + + + end if ! night or day + end associate + end subroutine LeafLayerPhotosynthesis + +! ======================================================================================= + function ft1_f(tl, ha) result(ans) ! !!DESCRIPTION: @@ -1157,4 +1245,146 @@ subroutine quadratic_f (a, b, c, r1, r2) end subroutine quadratic_f + ! ==================================================================================== + + subroutine UpdateCanopyNCanNRadPresent(currentPatch) + + ! --------------------------------------------------------------------------------- + ! This subroutine calculates two patch level quanities: + ! currentPatch%ncan and + ! currentPatch%present + ! + ! currentPatch%ncan(:,:) is a two dimensional array that indicates + ! the total number of leaf layers (including those that are not exposed to light) + ! in each canopy layer and for each functional type. + ! + ! currentPatch%nrad(:,:) is a two dimensional array that indicates + ! the total number of EXPOSED leaf layers, but for all intents and purposes + ! in the photosynthesis routine, this appears to be the same as %ncan... + ! + ! currentPatch%present(:,:) has the same dimensions, is binary, and + ! indicates whether or not leaf layers are present (by evaluating the canopy area + ! profile). + ! --------------------------------------------------------------------------------- + + use EDTypesMod, only : cp_nclmax + use EDTypesMOd, only : numpft_ed + + ! Arguments + type(ed_patch_type), target :: currentPatch + type(ed_cohort_type), pointer :: currentCohort + + ! Locals + integer :: CL ! Canopy Layer Index + integer :: ft ! Function Type Index + integer :: iv ! index of the exposed leaf layer for each canopy layer and pft + + ! Loop through the cohorts in this patch, associate each cohort with a layer and PFT + ! and use the cohort's memory of how many layer's it takes up to assign the maximum + ! of the layer/pft index it is in + ! --------------------------------------------------------------------------------- + + currentPatch%ncan(:,:) = 0 + !redo the canopy structure algorithm to get round a bug that is happening for site 125, FT13. + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + + currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & + max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) + + currentCohort => currentCohort%shorter + + enddo !cohort + + ! NRAD = NCAN ... + currentPatch%nrad = currentPatch%ncan + + ! Now loop through and identify which layer and pft combo has scattering elements + do CL = 1,cp_nclmax + do ft = 1,numpft_ed + currentPatch%present(CL,ft) = 0 + do iv = 1, currentPatch%nrad(CL,ft); + if(currentPatch%canopy_area_profile(CL,ft,iv) > 0._r8)then + currentPatch%present(CL,ft) = 1 + end if + end do !iv + enddo !ft + enddo !CL + + return + end subroutine UpdateCanopyNCanNRadPresent + + ! ==================================================================================== + + subroutine GetCanopyGasParameters(can_press, can_o2_partialpress, & + veg_temp,mm_kco2,mm_ko2,co2_comppoint) + + ! --------------------------------------------------------------------------------- + ! This subroutine calculates the specific Michaelis Menten Parameters (pa) for CO2 + ! and O2, as well as the CO2 compentation point. + ! --------------------------------------------------------------------------------- + + use FatesConstantsMod, only :: umol_per_mol + use FatesConstantsMod, only :: mmol_per_mol + + ! Arguments + real(r8), intent(in) :: can_press ! Air pressure within the canopy (Pa) + real(r8), intent(in) :: can_o2_partialpress ! Partial press of o2 in the canopy (Pa) + real(r8), intent(in) :: veg_tempk ! The temperature of the vegetation (K) + + real(r8), intent(out) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(out) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8), intent(out) :: co2_comppoint ! CO2 compensation point (Pa) + + ! Locals + real(r8) :: kc25 ! Michaelis-Menten constant for CO2 at 25C (Pa) + real(r8) :: ko25 ! Michaelis-Menten constant for O2 at 25C (Pa) + real(r8) :: sco ! relative specificity of rubisco + real(r8) :: cp25 ! CO2 compensation point at 25C (Pa) + + ! --------------------------------------------------------------------------------- + ! Intensive values (per mol of air) + ! kc, ko, currentPatch, from: Bernacchi et al (2001) + ! Plant, Cell and Environment 24:253-259 + ! --------------------------------------------------------------------------------- + + real(r8), parameter :: mm_kc25_umol_per_mol = 404.9_r8 + real(r8), parameter :: mm_ko25_mmol_per_mol = 278.4_r8 + real(r8), parameter :: co2_comppoint_umol_per_mol = 42.75_r8 + + ! Activation energy, from: + ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 + ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 + ! except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282 + + real(r8), parameter :: kcha = 79430._r8 ! activation energy for kc (J/mol) + real(r8), parameter :: koha = 36380._r8 ! activation energy for ko (J/mol) + real(r8), parameter :: cpha = 37830._r8 ! activation energy for cp (J/mol) + + + ! Derive sco from currentPatch and O2 using present-day O2 (0.209 mol/mol) and re-calculate + ! currentPatch to account for variation in O2 using currentPatch = 0.5 O2 / sco + + ! FIXME (RGK 11-30-2016 THere are more constants here, but I don't have enough information + ! about what they are or do, so I can't give them more descriptive names. Someone please + ! fill this in when possible) + + kc25 = ( mm_kc25_umol_per_mol / umol_per_mol ) * can_press + ko25 = ( mm_ko25_mmol_per_mol / mmol_per_mol ) * can_press + sco = 0.5_r8 * 0.209_r8 / (co2_comppoint_umol_per_mol / umol_per_mol ) + cp25 = 0.5_r8 * can_o2_partialpress / sco + + if( veg_tempk.gt.150_r8 .and. veg_tempk.lt.350_r8 )then + mm_kco2 = kc25 * ft1_f(veg_tempk, kcha) + mm_ko = ko25 * ft1_f(veg_tempk, koha) + co2_comppoint = cp25 * ft1_f(veg_tempk, cpha) + else + mm_kco2 = 1.0_r8 + mm_ko = 1.0_r8 + co2_comppoint = 1.0_r8 + end if + + return + end subroutine GetCanopyGasParameters + end module EDPhotosynthesisMod From 7258aceb42d6b1484453f985d65ca10360de34e6 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 1 Dec 2016 18:41:55 -0800 Subject: [PATCH 02/15] Photosynthesis refactor: most of the code has been re-organized. Next step is careful step-by-step pass through comparing with existing version. --- biogeophys/EDPhotosynthesisMod.F90 | 1163 +++++++++++++++------------- 1 file changed, 626 insertions(+), 537 deletions(-) diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index a0597d1ee5..f1dfe07ae6 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -1,30 +1,45 @@ module EDPhotosynthesisMod + + !------------------------------------------------------------------------------ + ! !DESCRIPTION: + ! Calculates the photosynthetic fluxes for the ED model + ! This code is equivalent to the 'photosynthesis' subroutine in PhotosynthesisMod.F90. + ! We have split this out to reduce merge conflicts until we can pull out + ! common code used in both the ED and CLM versions. + ! + ! Parameter for activation and deactivation energies were taken from: + ! Activation energy, from: + ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 + ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 + ! except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282 + ! High temperature deactivation, from: + ! Leuning (2002) Plant, Cell and Environment 25:1205-1210 + ! The factor "c" scales the deactivation to a value of 1.0 at 25C + ! ------------------------------------------------------------------------------------ + + ! !USES: + ! - !------------------------------------------------------------------------------ - ! !DESCRIPTION: - ! Calculates the photosynthetic fluxes for the ED model - ! This code is equivalent to the 'photosynthesis' subroutine in PhotosynthesisMod.F90. - ! We have split this out to reduce merge conflicts until we can pull out - ! common code used in both the ED and CLM versions. - ! - ! !USES: - ! - - 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: - public :: Photosynthesis_ED !ED specific photosynthesis routine - - character(len=*), parameter, private :: sourcefile = & - __FILE__ - !------------------------------------------------------------------------------ + 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: + public :: Photosynthesis_ED !ED specific photosynthesis routine + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + !------------------------------------------------------------------------------ + + ! maximum stomatal resistance [s/m] (used across several procedures) + real(r8),parameter :: rsmax0 = 2.e4_r8 + + logical :: DEBUG = .false. contains @@ -90,10 +105,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) integer , parameter :: psn_type = 2 !c3 or c4. - logical :: DEBUG = .false. - + real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) + ! either from cohort or patch-pft ! ! Leaf photosynthesis parameters + ! Note: None of these variables need to be an array. We put them + ! in arrays only to enable user debugging diagnostics real(r8) :: vcmax_z(cp_nclmax,mxpft,cp_nlevcan) ! maximum rate of carboxylation (umol co2/m**2/s) real(r8) :: jmax_z(cp_nclmax,mxpft,cp_nlevcan) ! maximum electron transport rate (umol electrons/m**2/s) real(r8) :: tpu_z(cp_nclmax,mxpft,cp_nlevcan) ! triose phosphate utilization rate (umol CO2/m**2/s) @@ -101,19 +118,19 @@ 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 ! intracellular leaf CO2 (Pa) - real(r8) :: lnc(mxpft) ! leaf N concentration (gN leaf/m^2) - - real(r8) :: kc ! Michaelis-Menten constant for CO2 (Pa) - real(r8) :: ko ! Michaelis-Menten constant for O2 (Pa) - real(r8) :: co2_cp ! CO2 compensation point (Pa) + real(r8) :: anet_av(cp_nclmax,mxpft,cp_nlevcan) ! net leaf photosynthesis (umol CO2/m**2/s) + ! averaged over sun and shade leaves. + + real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) + real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8) :: co2_cpoint ! CO2 compensation point (Pa) ! --------------------------------------------------------------- ! TO-DO: bbbopt is slated to be transferred to the parameter file ! ---------------------------------------------------------------- real(r8) :: bbbopt(psn_type) ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s) - real(r8) :: bbb(mxpft) ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) + real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation at 25C (umol CO2/m**2/s) @@ -122,83 +139,24 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8) :: lmr25top(mxpft) ! canopy top: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) real(r8) :: kp25top(mxpft) ! canopy top: initial slope of CO2 response curve (C4 plants) at 25C - real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C (umol CO2/m**2/s) - real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C (umol electrons/m**2/s) - real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) - real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - real(r8) :: kp25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C - - real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) - real(r8) :: jmaxha ! activation energy for jmax (J/mol) - real(r8) :: tpuha ! activation energy for tpu (J/mol) - real(r8) :: lmrha ! activation energy for lmr (J/mol) - - real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) - real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) - real(r8) :: tpuhd ! deactivation energy for tpu (J/mol) - real(r8) :: lmrhd ! deactivation energy for lmr (J/mol) - - real(r8) :: vcmaxse ! entropy term for vcmax (J/mol/K) - real(r8) :: jmaxse ! entropy term for jmax (J/mol/K) - real(r8) :: tpuse ! entropy term for tpu (J/mol/K) - real(r8) :: lmrse ! entropy term for lmr (J/mol/K) - - real(r8) :: vcmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8) :: jmaxc ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8) :: tpuc ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8) :: lmrc ! scaling factor for high temperature inhibition (25 C = 1.0) - - real(r8) :: qe(psn_type) ! quantum efficiency, used only for C4 (mol CO2 / mol photons) - real(r8) :: fnps ! fraction of light absorbed by non-photosynthetic pigments - real(r8) :: theta_psii ! empirical curvature parameter for electron transport rate - - real(r8) :: theta_cj(psn_type) ! empirical curvature parameter for ac, aj photosynthesis co-limitation - real(r8) :: theta_ip ! empirical curvature parameter for ap photosynthesis co-limitation - ! Other - integer :: c,CL,f,s,iv,j,ps,ft,ifp ! indices + integer :: cl,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) :: 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) - real(r8) :: gs_mol ! leaf stomatal conductance (umol H2O/m**2/s) - real(r8) :: gs ! leaf stomatal conductance (m/s) - real(r8) :: hs ! fractional humidity at leaf surface (dimensionless) - real(r8) :: tl ! leaf temperature in photosynthesis temperature function (K) - real(r8) :: ha ! activation energy in photosynthesis temperature function (J/mol) - real(r8) :: hd ! deactivation energy in photosynthesis temperature function (J/mol) - real(r8) :: se ! entropy term in photosynthesis temperature function (J/mol/K) - real(r8) :: cc2 ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8) :: ciold ! previous value of Ci for convergence check - real(r8) :: gs_mol_err ! gs_mol for error check - real(r8) :: je ! electron transport rate (umol electrons/m**2/s) - real(r8) :: qabs ! PAR absorbed by PS II (umol photons/m**2/s) - real(r8) :: aquad,bquad,cquad ! terms for quadratic equations - real(r8) :: r1,r2 ! roots of quadratic equation real(r8) :: ceair ! vapor pressure of air, constrained (Pa) - real(r8) :: act25 ! (umol/mgRubisco/min) Rubisco activity at 25 C - integer :: niter ! iteration loop index real(r8) :: nscaler ! leaf nitrogen scaling coefficient real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass +! real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) - real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) - real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) - real(r8) :: ap ! product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m**2/s) - - 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 + real(r8) :: laifrac real(r8) :: tcsoi ! Temperature response function for root respiration. - real(r8) :: tc ! Temperature response function for wood + real(r8) :: tcwood ! Temperature response function for wood - - real(r8) :: q10 ! temperature dependence of root respiration - integer :: sunsha ! sun (1) or shaded (2) leaves... - real(r8) :: coarse_wood_frac ! amount of woody biomass that is coarse... +! real(r8) :: coarse_wood_frac ! amount of woody biomass that is coarse... real(r8) :: tree_area real(r8) :: gs_cohort real(r8) :: rscanopy @@ -222,84 +180,25 @@ 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_r8 - real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8 - - - associate( & - c3psn => pftcon%c3psn , & ! photosynthetic pathway: 0. = c4, 1. = c3 + associate( & + c3psn => pftcon%c3psn , & slatop => pftcon%slatop , & ! specific leaf area at top of canopy, projected area basis [m^2/gC] flnr => pftcon%flnr , & ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) 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 + frootcn => pftcon%frootcn , & ! froot C:N (gc/gN) ! slope of BB relationship + q10 => EDParamsShareInst%Q10) - ! 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 - ! simulatoins - q10 = 1.5_r8 - Q10 = EDParamsShareInst%Q10 !==============================================================================! ! Photosynthesis and stomatal conductance parameters, from: ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 !==============================================================================! - ! vcmax25 parameters, from CN - - act25 = 3.6_r8 !umol/mgRubisco/min - ! Convert rubisco activity units from umol/mgRubisco/min -> umol/gRubisco/s - act25 = act25 * mg_per_g / sec_per_min - - ! Activation energy, from: - ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 - ! Bernacchi et al (2003) Plant, Cell and Environment 26:1419-1430 - ! except TPU from: Harley et al (1992) Plant, Cell and Environment 15:271-282 - - vcmaxha = 65330._r8 - jmaxha = 43540._r8 - tpuha = 53100._r8 - lmrha = 46390._r8 - - ! High temperature deactivation, from: - ! Leuning (2002) Plant, Cell and Environment 25:1205-1210 - ! The factor "c" scales the deactivation to a value of 1.0 at 25C - - vcmaxhd = 149250._r8 - jmaxhd = 152040._r8 - tpuhd = 150650._r8 - lmrhd = 150650._r8 - - vcmaxse = 485._r8 - jmaxse = 495._r8 - tpuse = 490._r8 - lmrse = 490._r8 - - vcmaxc = fth25_f(vcmaxhd, vcmaxse) - jmaxc = fth25_f(jmaxhd, jmaxse) - tpuc = fth25_f(tpuhd, tpuse) - lmrc = fth25_f(lmrhd, lmrse) - ! Miscellaneous parameters, from Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - fnps = 0.15_r8 - theta_psii = 0.7_r8 - theta_ip = 0.95_r8 - - qe(1) = 0._r8 - theta_cj(1) = 0.98_r8 bbbopt(1) = 10000._r8 - - qe(2) = 0.05_r8 - theta_cj(2) = 0.80_r8 bbbopt(2) = 40000._r8 do s = 1,nsites @@ -341,24 +240,26 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Part IV. Identify some environmentally derived parameters: + ! These quantities are biologically irrelevant ! Michaelis-Menten constant for CO2 (Pa) ! Michaelis-Menten constant for O2 (Pa) ! CO2 compensation point (Pa) - - call GetCanopyGasParameters(bc_in(s)%forc_pbot,bc_in(s)%oair_pa(ifp), & - bc_in(s)%t_veg_pa(ifp),kc,ko,co2_cp) - - - ! 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 - - ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures - ! that hs does not go to zero. Also eair <= esat_tv so that hs <= 1 - ceair = min( max(bc_in(s)%eair_pa(ifp), 0.05_r8*bc_in(s)%esat_tv_pa(ifp)), bc_in(s)%esat_tv_pa(ifp) ) - - + ! CF? I have no idea what cf is (rgk 12-01-2016) + ! leaf boundary layer conductance of h20 + ! constrained vapor pressure + call GetCanopyGasParameters(bc_in(s)%forc_pbot, & ! in + bc_in(s)%oair_pa(ifp), & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%tgcm_pa(ifp), & ! in + bc_in(s)%eair_pa(ifp), & ! in + bc_in(s)%esat_tv_pa(ifp), & ! in + bc_in(s)%rb_pa(ifp), & ! in + mm_kco2, & ! out + mm_ko2, & ! out + co2_cpoint, & ! out + cf, & ! out + gb_mol, & ! out + ceair) ! out ! Part V. Pre-process some variables that are PFT dependent ! but not environmentally dependent @@ -367,7 +268,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) do FT = 1,numpft_ed !calculate patch and pft specific properties at canopy top. ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) - lnc(FT) = 1._r8 / (slatop(FT) * leafcn(FT)) + lnc = 1._r8 / (slatop(FT) * leafcn(FT)) !at the moment in ED we assume that there is no active N cycle. This should change, of course. FIX(RF,032414) Sep2011. vcmax25top(FT) = fnitr(FT) !fudge - shortcut using fnitr as a proxy for vcmax... @@ -403,7 +304,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) / (umolC_to_kgC * g_per_kg) + lmr25top(FT) = lmr25top(FT) * lnc / (umolC_to_kgC * g_per_kg) end do !FT @@ -412,7 +313,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! variables will be available at the cohort scale, and not the ! pft scale. So here we split and use different looping structures ! ------------------------------------------------------------------ -! if ( use_fates_plant_hydro ) + ! if ( use_fates_plant_hydro ) !==============================================================================! @@ -433,13 +334,13 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! calculated every timestep. Others are calculated only if daytime do iv = 1, currentPatch%nrad(CL,FT) - if (use_fates_plant_hydro) then - !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) - !! btran = currentCohort%btran(iv) - else - !! bbb = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8) - !! btran = currentPatch%btran_ft(currentCohort%pft) - end if + !! if (use_fates_plant_hydro) then + !! !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) + !! !! btran = currentCohort%btran(iv) + !! else + bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(FT), 1._r8) + btran_eff = currentPatch%btran_ft(ft) + !! end if vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) !vegetation area index. if (iv == 1) then @@ -451,42 +352,64 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Scale for leaf nitrogen profile nscaler = exp(-kn(FT) * laican) - - - call LeafLayerPhotosynthesis(currentPatch%ed_parsun_z(CL,FT,iv), & ! in - currentPatch%ed_parsha_z(CL,FT,iv), & ! in - currentPatch%ed_laisun_z(CL,FT,iv), & ! in - currentPatch%ed_laisha_z(CL,FT,iv), & ! in - currentPatch%canopy_area_profile(CL,FT,iv), & ! in - ft, & ! in - nscaler, & ! in - lmr25top(ft), & ! in - vcmax25top(ft), & ! in - jmax25top(ft), & ! in - tpu25top(ft), & ! in - kp25top(ft), & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - btran, & ! in - bbb, & ! in - cf, & ! in - gb_mol, & ! in - lmr_z(CL,FT,iv), & ! out - currentPatch%psn_z(cl,ft,iv), & ! out - rs_z(CL,FT,iv), & ! out - an_av(CL,FT,iv)) ! out - - + + call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in + nscaler, & ! in + ft, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + lmr_z(CL,FT,iv)) ! out + call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(CL,FT,iv), & ! in + ft, & ! in + vcmax25top(ft), & ! in + jmax25top(ft), & ! in + tpu25top(ft), & ! in + kp25top(ft), & ! in + nscaler, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + btran_eff, & ! in + vcmax_z(CL,FT,iv), & ! out + jmax_z(CL,FT,iv), & ! out + tpu_z(CL,FT,iv), & ! out + kp_z(CL,FT,iv) ) ! out + + call LeafLayerPhotosynthesis(currentPatch%f_sun(CL,FT,iv), & ! in + currentPatch%ed_parsun_z(CL,FT,iv), & ! in + currentPatch%ed_parsha_z(CL,FT,iv), & ! in + currentPatch%ed_laisun_z(CL,FT,iv), & ! in + currentPatch%ed_laisha_z(CL,FT,iv), & ! in + currentPatch%canopy_area_profile(CL,FT,iv), & ! in + ft, & ! in + nscaler, & ! in + vcmax_z(CL,FT,iv), & ! in + jmax_z(CL,FT,iv), & ! in + tpu_z(CL,FT,iv), & ! in + kp_z(CL,FT,iv), & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%esat_tv_pa(ifp), & ! in + bc_in(s)%forc_pbot, & ! in + bc_in(s)%cair_pa(ifp), & ! in + bc_in(s)%oair_pa(ifp), & ! in + btran_eff, & ! in + bbb, & ! in + cf, & ! in + gb_mol, & ! in + ceair, & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr_z(CL,FT,iv), & ! in + currentPatch%psn_z(cl,ft,iv), & ! out + rs_z(CL,FT,iv), & ! out + anet_av(CL,FT,iv)) ! out + end do ! iv end if !present enddo !PFT enddo !CL - ! Leaf boundary layer conductance, umol/m**2/s - - !==============================================================================! ! Unpack fluxes from arrays into cohorts @@ -526,14 +449,14 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) if (currentCohort%nv > 1) then !is there canopy, and are the leaves on? - currentCohort%gpp_tstep = sum(currentPatch%psn_z(cl,ft,1:currentCohort%nv-1) * & + 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%rdark = 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) * umolC_to_kgC * dtime + 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) = anet_av(cl,ft,1:currentCohort%nv) * umolC_to_kgC * dtime else @@ -611,9 +534,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! 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) + tcwood = 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_stem_n * base_mr_20 * tc + currentCohort%livestem_mr = live_stem_n * base_mr_20 * tcwood else currentCohort%livestem_mr = 0._r8 end if @@ -728,28 +651,36 @@ end subroutine Photosynthesis_ED ! ======================================================================================= -subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in +subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in + parsun_lsl, & ! in parsha_lsl, & ! in laisun_lsl, & ! in laisha_lsl, & ! in canopy_area_lsl, & ! in ft, & ! in nscaler, & ! in - lmr25top_ft, & ! in - vcmax25top_ft, & ! in - jmax25top_ft, & ! in - tpu25top_ft, & ! in - kp25top_ft, & ! in - t_veg, & ! in + vcmax, & ! in + jmax, & ! in + tpu, & ! in + co2_rcurve_islope, & ! in + veg_tempk, & ! in + veg_esat, & ! in + can_press, & ! in + can_co2_ppress, & ! in + can_o2_ppress, & ! in btran, & ! in bbb, & ! in cf, & ! in gb_mol, & ! in - lmr_out, & ! out + ceair, & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr, & ! in psn_out, & ! out rstoma_out, & ! out anet_av_out) ! out - + ! ------------------------------------------------------------------------------------ ! This subroutine calculates photosynthesis and stomatal conductance within each leaf ! sublayer. @@ -759,141 +690,111 @@ subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in ! Other arguments or variables may be indicative of scales broader than the LSL. ! ------------------------------------------------------------------------------------ + use EDEcophysContype , only : EDecophyscon + use pftconMod , only : pftcon + ! Arguments ! ------------------------------------------------------------------------ + real(r8), intent(in) :: f_sun_lsl real(r8), intent(in) :: parsun_lsl real(r8), intent(in) :: parsha_lsl real(r8), intent(in) :: laisun_lsl real(r8), intent(in) :: laisha_lsl - real(r8), intent(in) :: elai_lsl - real(r8), intent(in) :: esai_lsl real(r8), intent(in) :: canopy_area_lsl - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: lmr25top_ft ! canopy top leaf maint resp rate at 25C - ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C - ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C - ! for this pft (umol electrons/m**2/s) - real(r8), intent(in) :: tpu25top_ft ! canopy top triose phosphate utilization rate at 25C - ! for this pft (umol CO2/m**2/s) - real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve - ! (C4 plants) at 25C, canopy top, this pft - real(r8), intent(in) :: t_veg ! vegetation temperature - real(r8), intent(in) :: btran ! - real(r8), intent(in) :: bbb - real(r8), intent(in) :: cf - real(r8), intent(in) :: gb_mol + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(in) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) + real(r8), intent(in) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) + + ! Important Note on the following gas pressures. This photosynthesis scheme will iteratively + ! solve for the co2 partial pressure at the leaf surface (ie in the stomata). The reference + ! point for these input values are NOT within that boundary layer that separates the stomata from + ! the canopy air space. The reference point for these is on the outside of that boundary + ! layer. This routine, which operates at the leaf scale, makes no assumptions about what the + ! scale of the refernce is, it could be lower atmosphere, it could be within the canopy + ! but most likely it is the closest value one can get to the edge of the leaf's boundary + ! layer. We use the convention "can_" because a reference point of within the canopy + ! ia a best reasonable scenario of where we can get that information from. + + real(r8), intent(in) :: can_press ! Air pressure NEAR the surface of the leaf (Pa) + real(r8), intent(in) :: can_co2_ppress ! Partial pressure of CO2 NEAR the leaf surface (Pa) + real(r8), intent(in) :: can_o2_ppress ! Partial pressure of O2 NEAR the leaf surface (Pa) + real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) + real(r8), intent(in) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) + real(r8), intent(in) :: cf ! s m**2/umol -> s/m + real(r8), intent(in) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8), intent(in) :: ceair ! vapor pressure of air, constrained (Pa) + real(r8), intent(in) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8), intent(in) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8), intent(in) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8), intent(in) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) - real(r8), intent(out) :: lmr_out - real(r8), intent(out) :: psn_out - real(r8), intent(out) :: rstoma_out ! stomatal resistance (1/gs_lsl) (s/m) - real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) averaged over sun and shade leaves. - real(r8), intent(out) :: gstoma_out ! Stomatal Conductance of this leaf layer (m/s) - + real(r8), intent(out) :: psn_out ! carbon assimilated in this leaf layer umolC/m2/s + real(r8), intent(out) :: rstoma_out ! stomatal resistance (1/gs_lsl) (s/m) + real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) + ! averaged over sun and shade leaves. ! Locals ! ------------------------------------------------------------------------ - integer :: ps ! Index for the different photosynthetic pathways C3,C4 - integer :: sunsha ! Index for differentiating sun and shade - - real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C (umol CO2/m**2/s) - real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C (umol electrons/m**2/s) - real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) - real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C - - - !!!!!!! NOTE: WE CAN REMOVE THE LSL IDENTIFIER FROM LOCALS AFTER THEY HAVE MORE VERBOSE NAMES - ! - - real(r8) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) - real(r8) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) - - + integer :: ps ! Index for the different photosynthetic pathways C3,C4 + integer :: sunsha ! Index for differentiating sun and shade + real(r8) :: gstoma ! Stomatal Conductance of this leaf layer (m/s) real(r8) :: agross ! co-limited gross leaf photosynthesis (umol CO2/m**2/s) real(r8) :: anet ! net leaf photosynthesis (umol CO2/m**2/s) + real(r8) :: je ! electron transport rate (umol electrons/m**2/s) + real(r8) :: qabs ! PAR absorbed by PS II (umol photons/m**2/s) + real(r8) :: aquad,bquad,cquad ! terms for quadratic equations + real(r8) :: r1,r2 ! roots of quadratic equation + real(r8) :: co2_intra_c ! intracellular leaf CO2 (Pa) + real(r8) :: co2_intra_c_old ! intracellular leaf CO2 (Pa) (previous iteration) + logical :: loop_continue ! Loop control variable + integer :: niter ! iteration loop index + real(r8) :: gs_mol ! leaf stomatal conductance (umol H2O/m**2/s) + real(r8) :: gs ! leaf stomatal conductance (m/s) + real(r8) :: hs ! fractional humidity at leaf surface (dimensionless) + real(r8) :: gs_mol_err ! gs_mol for error check + real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) + real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) + real(r8) :: ap ! product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m**2/s) + real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) + real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + ! Parameters + ! ------------------------------------------------------------------------ + ! Fraction of light absorbed by non-photosynthetic pigments + real(r8),parameter :: fnps = 0.15_r8 + + ! empirical curvature parameter for electron transport rate + real(r8),parameter :: theta_psii = 0.7_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_r8 + real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8 + ! quantum efficiency, used only for C4 (mol CO2 / mol photons) + real(r8),parameter,dimension(2) :: quant_eff = [0.0_r8,0.05_r8] + + ! empirical curvature parameter for ac, aj photosynthesis co-limitation + real(r8),parameter,dimension(2) :: theta_cj = [0.98_r8,0.80_r8] + + ! empirical curvature parameter for ap photosynthesis co-limitation + real(r8),parameter :: theta_ip = 0.95_r8 - ! Parameters - ! ------------------------------------------------------------------------ - associate( & - c3psn => pftcon%c3psn , & ! photosynthetic pathway: 0. = c4, 1. = c3 - slatop => pftcon%slatop , & ! specific leaf area at top of canopy, projected area basis [m^2/gC] - flnr => pftcon%flnr , & ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) - 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 - - + associate( c3psn => pftcon%c3psn, & ! photosynthetic pathway: 0. = c4, 1. = c3 + bb_slope => EDecophyscon%BB_slope ) ! slope of BB relationship + if (nint(c3psn(ft)) == 1)then ps = 1 else ps = 2 end if - - - ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s - ! ---------------------------------------------------------------------------------- - lmr25 = lmr25top_ft * nscaler - if ( nint(c3psn(ft)) == 1)then - lmr_out = lmr25 * ft1_f(t_veg, lmrha) * & - fth_f(t_veg, lmrhd, lmrse, lmrc) - else - lmr_out = lmr25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) - lmr_out = lmr_out / (1._r8 + exp( 1.3_r8*(t_veg-(tfrz+55._r8)) )) - end if - - - - - ! Part II: Localized Biophysical Rates - ! ---------------------------------------------------------------------------------- - - if ( parsun_lsl <= 0._r8) then ! night time - vcmax = 0._r8 - jmax = 0._r8 - tpu = 0._r8 - co2_rcurve_islope = 0._r8 - else ! day time - vcmax25 = vcmax25top_ft * nscaler - jmax25 = jmax25top_ft * nscaler - tpu25 = tpu25top_ft * nscaler - co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler - - ! Adjust for temperature - vcmax = vcmax25 * ft1_f(t_veg, vcmaxha) * fth_f(t_veg, vcmaxhd, vcmaxse, vcmaxc) - jmax = jmax25 * ft1_f(t_veg, jmaxha) * fth_f(t_veg, jmaxhd, jmaxse, jmaxc) - tpu = tpu25 * ft1_f(t_veg, tpuha) * fth_f(t_veg, tpuhd, tpuse, tpuc) - - if (nint(c3psn(FT)) /= 1) then - vcmax = vcmax25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) - vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-t_veg ) )) - vcmax = vcmax / (1._r8 + exp( 0.3_r8*(t_veg-(tfrz+40._r8)) )) - end if - co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((t_veg-(tfrz+25._r8))/10._r8) !q10 response of product limited psn. - end if - - ! Adjust for water limitations - vcmax = vcmax * btran - - ! Leaf Maintenance Respiration has no direct water limitation effect - ! lmr_out = lmr_out * (nothing) - - - - - - ! Part III: Photosynthesis and Conductance ! ---------------------------------------------------------------------------------- @@ -901,7 +802,7 @@ subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in anet_av_out = 0._r8 psn_out = 0._r8 - rstoma_out = min(rsmax0, 1._r8/bbb * cf) + rstoma_out = min(rsmax0, 1._r8/bbb * cf) else ! day time (a little bit more complicated ...) @@ -915,9 +816,9 @@ subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. - rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. + rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. anet_av_out = 0._r8 - gstoma_out = 0._r8 + gstoma = 0._r8 do sunsha = 1,2 ! Electron transport rate for C3 plants. Convert par from W/m2 to umol photons/m**2/s @@ -927,7 +828,7 @@ subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in if(sunsha == 1)then !sunlit if(( laisun_lsl * canopy_area_lsl) > 0.0000000001_r8)then - qabs = parsun_lsl / (laisun_lsl * canopy_area_profile_lsl ) + qabs = parsun_lsl / (laisun_lsl * canopy_area_lsl ) qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 else @@ -935,182 +836,183 @@ subroutine LeafLayerPhotosynthesis(parsun_lsl, & ! in end if else - qabs = currentPatch%ed_parsha_lsl / (currentPatch%ed_laisha_lsl * & - currentPatch%canopy_area_profile_lsl) - qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 - - end if - - !convert the absorbed par into absorbed par per m2 of leaf, - ! so it is consistant with the vcmax and lmr numbers. - aquad = theta_psii - bquad = -(qabs + jmax_z(cl,ft,iv)) - cquad = qabs * jmax_z(cl,ft,iv) - call quadratic_f (aquad, bquad, cquad, r1, r2) - je = min(r1,r2) - - ! Iterative loop for ci beginning with initial guess - ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) - - if (nint(c3psn(FT)) == 1)then - ci = init_a2l_co2_c3 * bc_in(s)%cair_pa(ifp) - else - ci = init_a2l_co2_c4 * bc_in(s)%cair_pa(ifp) - end if - - niter = 0 - exitloop = 0 - do while(exitloop == 0) - ! Increment iteration counter. Stop if too many iterations - niter = niter + 1 - - ! Save old ci - ciold = ci - - ! Photosynthesis limitation rate calculations - if (nint(c3psn(FT)) == 1)then - ! C3: Rubisco-limited photosynthesis - ac = vcmax_z(cl,ft,iv) * max(ci-co2_cp, 0._r8) / (ci+kc * & - (1._r8+bc_in(s)%oair_pa(ifp)/ko)) - ! C3: RuBP-limited photosynthesis - aj = je * max(ci-co2_cp, 0._r8) / (4._r8*ci+8._r8*co2_cp) - ! C3: Product-limited photosynthesis - ap = 3._r8 * tpu_z(cl,ft,iv) - else - ! C4: Rubisco-limited photosynthesis - ac = vcmax_z(cl,ft,iv) - ! C4: RuBP-limited photosynthesis - if(sunsha == 1)then !sunlit - if((currentPatch%ed_laisun_z(cl,ft,iv) * currentPatch%canopy_area_profile(cl,ft,iv)) > & - 0.0000000001_r8)then !guard against /0's in the night. - aj = qe(ps) * currentPatch%ed_parsun_z(cl,ft,iv) * 4.6_r8 - !convert from per cohort to per m2 of leaf) - aj = aj / (currentPatch%ed_laisun_z(cl,ft,iv) * & - currentPatch%canopy_area_profile(cl,ft,iv)) - else - aj = 0._r8 - end if - else - aj = qe(ps) * currentPatch%ed_parsha_z(cl,ft,iv) * 4.6_r8 - aj = aj / (currentPatch%ed_laisha_z(cl,ft,iv) * & - currentPatch%canopy_area_profile(cl,ft,iv)) - end if - - ! C4: PEP carboxylase-limited (CO2-limited) - ap = co2_rcurve_islope_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) - bquad = -(ac + aj) - cquad = ac * aj - call quadratic_f (aquad, bquad, cquad, r1, r2) - ai = min(r1,r2) - - aquad = theta_ip - bquad = -(ai + ap) - cquad = ai * ap - call quadratic_f (aquad, bquad, cquad, r1, r2) - agross = min(r1,r2) - - ! Net carbon assimilation. Exit iteration if an < 0 - anet = agross - lmr_out - if (anet < 0._r8) then - exitloop = 1 - end if - - ! Quadratic gs_mol calculation with an known. Valid for an >= 0. - ! With an <= 0, then gs_mol = bbb - - cs = bc_in(s)%cair_pa(ifp) - 1.4_r8/gb_mol * anet * bc_in(s)%forc_pbot - cs = max(cs,1.e-06_r8) - aquad = cs - bquad = cs*(gb_mol - bbb(FT)) - bb_slope(ft)*anet*bc_in(s)%forc_pbot - cquad = -gb_mol*(cs*bbb(FT) + & - bb_slope(ft)*anet*bc_in(s)%forc_pbot*ceair/bc_in(s)%esat_tv_pa(ifp)) - call quadratic_f (aquad, bquad, cquad, r1, r2) - gs_mol = max(r1,r2) - - ! Derive new estimate for ci - ci = bc_in(s)%cair_pa(ifp) - anet * 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-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 - - ! End of ci iteration. Check for an < 0, in which case gs_mol = bbb - if (anet < 0._r8) then - gs_mol = bbb(FT) - end if - - ! 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 * anet * bc_in(s)%forc_pbot - cs = max(cs,1.e-06_r8) - ci = bc_in(s)%cair_pa(ifp) - & - anet * 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 - - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 737 ', currentPatch%psn_z(cl,ft,iv) - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross - 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 - - psn_out = psn_out + agross * f_sun_lsl - anet_av_out = anet_av_out + anet * f_sun_lsl - gstoma_out = gstoma_out + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl - - else - - psn_out = psn_out + agross * (1.0_r8-f_sun_lsl) - anet_av_out = anet_av + anet * (1.0_r8-f_sun_lsl) - gstoma_out = gstoma_out + & - 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) - - end if - - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', currentPatch%psn_z(cl,ft,iv) - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', agross - 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 (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 - hs = (gb_mol*ceair + gs_mol*bc_in(s)%esat_tv_pa(ifp)) / ((gb_mol+gs_mol)*bc_in(s)%esat_tv_pa(ifp)) - gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/cs*bc_in(s)%forc_pbot + bbb(FT) - - if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then - write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' - write (fates_log(),*) gs_mol, gs_mol_err - end if - - enddo !sunsha loop - - !average leaf-level stomatal resistance rate over sun and shade leaves... - rstoma_out = 1._r8/gstoma_out - - else !No leaf area. This layer is present only because of stems. (leaves are off, or have reduced to 0 - currentPatch%psn_z(cl,ft,iv) = 0._r8 - rs_lsl = min(rsmax0, 1._r8/bbb(FT) * cf) - - end if !is there leaf area? - - - end if ! night or day - end associate - end subroutine LeafLayerPhotosynthesis + qabs = parsha_lsl / (laisha_lsl * canopy_area_lsl) + qabs = qabs * 0.5_r8 * (1._r8 - fnps) * 4.6_r8 + + end if + + !convert the absorbed par into absorbed par per m2 of leaf, + ! so it is consistant with the vcmax and lmr numbers. + aquad = theta_psii + bquad = -(qabs + jmax) + cquad = qabs * jmax + call quadratic_f (aquad, bquad, cquad, r1, r2) + je = min(r1,r2) + + ! Iterative loop for ci beginning with initial guess + ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) + + if (nint(c3psn(FT)) == 1)then + co2_intra_c = init_a2l_co2_c3 * can_co2_ppress + else + co2_intra_c = init_a2l_co2_c4 * can_co2_ppress + end if + + niter = 0 + loop_continue = .true. + do while(loop_continue) + ! Increment iteration counter. Stop if too many iterations + niter = niter + 1 + + ! Save old co2_intra_c + co2_intra_c_old = co2_intra_c + + ! Photosynthesis limitation rate calculations + if (nint(c3psn(FT)) == 1)then + + ! C3: Rubisco-limited photosynthesis + ac = vcmax * max(co2_intra_c-co2_cpoint, 0._r8) / & + (co2_intra_c+mm_kco2 * (1._r8+can_o2_ppress / mm_ko2 )) + + ! C3: RuBP-limited photosynthesis + aj = je * max(co2_intra_c-co2_cpoint, 0._r8) / (4._r8*co2_intra_c+8._r8*co2_cpoint) + + ! C3: Product-limited photosynthesis + ap = 3._r8 * tpu + + else + + ! C4: Rubisco-limited photosynthesis + ac = vcmax + + ! C4: RuBP-limited photosynthesis + if(sunsha == 1)then !sunlit + if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then !guard against /0's in the night. + aj = quant_eff(ps) * parsun_lsl * 4.6_r8 + !convert from per cohort to per m2 of leaf) + aj = aj / (laisun_lsl * canopy_area_lsl) + else + aj = 0._r8 + end if + else + aj = quant_eff(ps) * parsha_lsl * 4.6_r8 + aj = aj / (laisha_lsl * canopy_area_lsl) + end if + + ! C4: PEP carboxylase-limited (CO2-limited) + ap = co2_rcurve_islope * max(co2_intra_c, 0._r8) / can_press + + end if + + ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap + aquad = theta_cj(ps) + bquad = -(ac + aj) + cquad = ac * aj + call quadratic_f (aquad, bquad, cquad, r1, r2) + ai = min(r1,r2) + + aquad = theta_ip + bquad = -(ai + ap) + cquad = ai * ap + call quadratic_f (aquad, bquad, cquad, r1, r2) + agross = min(r1,r2) + + ! Net carbon assimilation. Exit iteration if an < 0 + anet = agross - lmr + if (anet < 0._r8) then + loop_continue = .false. + end if + + ! Quadratic gs_mol calculation with an known. Valid for an >= 0. + ! With an <= 0, then gs_mol = bbb + + leaf_co2_ppress = can_co2_ppress- 1.4_r8/gb_mol * anet * can_press + leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) + aquad = leaf_co2_ppress + bquad = leaf_co2_ppress*(gb_mol - bbb) - bb_slope(ft) * anet * can_press + cquad = -gb_mol*(leaf_co2_ppress*bbb + bb_slope(ft)*anet*can_press * ceair/ veg_esat ) + + call quadratic_f (aquad, bquad, cquad, r1, r2) + gs_mol = max(r1,r2) + + ! Derive new estimate for co2_intra_c + co2_intra_c = can_co2_ppress - anet * can_press * & + (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) + + ! Check for co2_intra_c convergence. Delta co2_intra_c/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(co2_intra_c-co2_intra_c_old)/can_press*1.e06_r8 <= 2.e-06_r8) .or. niter == 5) then + loop_continue = .false. + end if + end do !iteration loop + + ! End of co2_intra_c iteration. Check for an < 0, in which case gs_mol = bbb + if (anet < 0._r8) then + gs_mol = bbb + end if + + ! Final estimates for leaf_co2_ppress and co2_intra_c (needed for early exit of co2_intra_c iteration when an < 0) + leaf_co2_ppress = can_co2_ppress - 1.4_r8/gb_mol * anet * can_press + leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) + co2_intra_c = can_co2_ppress - anet * can_press * (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 + + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 737 ', psn_out + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl + + !accumulate total photosynthesis umol/m2 ground/s-1. weight per unit sun and sha leaves. + if(sunsha == 1)then !sunlit + psn_out = psn_out + agross * f_sun_lsl + anet_av_out = anet_av_out + anet * f_sun_lsl + gstoma = gstoma + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl + else + psn_out = psn_out + agross * (1.0_r8-f_sun_lsl) + anet_av_out = anet_av_out + anet * (1.0_r8-f_sun_lsl) + gstoma = gstoma + & + 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) + end if + + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', psn_out + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', agross + if ( DEBUG ) write(fates_log(),*) 'EDPhoto 760 ', f_sun_lsl + + ! Make sure iterative solution is correct + if (gs_mol < 0._r8) then + write (fates_log(),*)'Negative stomatal conductance:' + write (fates_log(),*)'gs_mol= ',gs_mol + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Compare with Ball-Berry model: gs_mol = m * an * hs/leaf_co2_ppress p + b + hs = (gb_mol*ceair + gs_mol* veg_esat ) / ((gb_mol+gs_mol)*veg_esat ) + gs_mol_err = bb_slope(ft)*max(anet, 0._r8)*hs/leaf_co2_ppress*can_press + bbb + + if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then + write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' + write (fates_log(),*) gs_mol, gs_mol_err + end if + + enddo !sunsha loop + + !average leaf-level stomatal resistance rate over sun and shade leaves... + rstoma_out = 1._r8/gstoma + + else !No leaf area. This layer is present only because of stems. (leaves are off, or have reduced to 0 + psn_out = 0._r8 + rstoma_out = min(rsmax0, 1._r8/bbb * cf) + + end if !is there leaf area? + + + end if ! night or day + end associate + return +end subroutine LeafLayerPhotosynthesis ! ======================================================================================= @@ -1267,9 +1169,11 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! profile). ! --------------------------------------------------------------------------------- - use EDTypesMod, only : cp_nclmax - use EDTypesMOd, only : numpft_ed - + use EDTypesMod , only : cp_nclmax + use EDTypesMOd , only : numpft_ed + use EDTypesMod , only : ed_patch_type + use EDTypesMod , only : ed_cohort_type + ! Arguments type(ed_patch_type), target :: currentPatch type(ed_cohort_type), pointer :: currentCohort @@ -1316,25 +1220,44 @@ end subroutine UpdateCanopyNCanNRadPresent ! ==================================================================================== - subroutine GetCanopyGasParameters(can_press, can_o2_partialpress, & - veg_temp,mm_kco2,mm_ko2,co2_comppoint) + subroutine GetCanopyGasParameters(can_press, & + can_o2_partialpress, & + veg_tempk, & + air_tempk, & + air_vpress, & + veg_esat, & + rb, & + mm_kco2, & + mm_ko2, & + co2_cpoint, & + cf, & + gb_mol, & + ceair) ! --------------------------------------------------------------------------------- ! This subroutine calculates the specific Michaelis Menten Parameters (pa) for CO2 ! and O2, as well as the CO2 compentation point. ! --------------------------------------------------------------------------------- - use FatesConstantsMod, only :: umol_per_mol - use FatesConstantsMod, only :: mmol_per_mol - + use FatesConstantsMod, only: umol_per_mol + use FatesConstantsMod, only: mmol_per_mol + use FatesConstantsMod, only : rgas => rgas_J_K_kmol + ! Arguments real(r8), intent(in) :: can_press ! Air pressure within the canopy (Pa) real(r8), intent(in) :: can_o2_partialpress ! Partial press of o2 in the canopy (Pa) real(r8), intent(in) :: veg_tempk ! The temperature of the vegetation (K) - + real(r8), intent(in) :: air_tempk ! Temperature of canopy air (K) + real(r8), intent(in) :: air_vpress ! Vapor pressure of canopy air (Pa) + real(r8), intent(in) :: veg_esat ! Saturated vapor pressure at veg surf (Pa) + real(r8), intent(in) :: rb ! Leaf Boundary layer resistance (s/m) + real(r8), intent(out) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) real(r8), intent(out) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) - real(r8), intent(out) :: co2_comppoint ! CO2 compensation point (Pa) + real(r8), intent(out) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8), intent(out) :: cf ! s m**2/umol -> s/m + real(r8), intent(out) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8), intent(out) :: ceair ! vapor pressure of air, constrained (Pa) ! Locals real(r8) :: kc25 ! Michaelis-Menten constant for CO2 at 25C (Pa) @@ -1350,7 +1273,7 @@ subroutine GetCanopyGasParameters(can_press, can_o2_partialpress, & real(r8), parameter :: mm_kc25_umol_per_mol = 404.9_r8 real(r8), parameter :: mm_ko25_mmol_per_mol = 278.4_r8 - real(r8), parameter :: co2_comppoint_umol_per_mol = 42.75_r8 + real(r8), parameter :: co2_cpoint_umol_per_mol = 42.75_r8 ! Activation energy, from: ! Bernacchi et al (2001) Plant, Cell and Environment 24:253-259 @@ -1371,20 +1294,186 @@ subroutine GetCanopyGasParameters(can_press, can_o2_partialpress, & kc25 = ( mm_kc25_umol_per_mol / umol_per_mol ) * can_press ko25 = ( mm_ko25_mmol_per_mol / mmol_per_mol ) * can_press - sco = 0.5_r8 * 0.209_r8 / (co2_comppoint_umol_per_mol / umol_per_mol ) + sco = 0.5_r8 * 0.209_r8 / (co2_cpoint_umol_per_mol / umol_per_mol ) cp25 = 0.5_r8 * can_o2_partialpress / sco if( veg_tempk.gt.150_r8 .and. veg_tempk.lt.350_r8 )then mm_kco2 = kc25 * ft1_f(veg_tempk, kcha) - mm_ko = ko25 * ft1_f(veg_tempk, koha) - co2_comppoint = cp25 * ft1_f(veg_tempk, cpha) + mm_ko2 = ko25 * ft1_f(veg_tempk, koha) + co2_cpoint = cp25 * ft1_f(veg_tempk, cpha) else mm_kco2 = 1.0_r8 - mm_ko = 1.0_r8 - co2_comppoint = 1.0_r8 + mm_ko2 = 1.0_r8 + co2_cpoint = 1.0_r8 end if + ! THESE HARD CODED CONVERSIONS NEED TO BE CALLED FROM GLOBAL CONSTANTS (RGK 10-13-2016) + cf = can_press/(rgas*1.e-3_r8 * air_tempk )*1.e06_r8 + gb_mol = (1._r8/ rb) * cf + + ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures + ! that hs does not go to zero. Also eair <= veg_esat so that hs <= 1 + ceair = min( max(air_vpress, 0.05_r8*veg_esat ),veg_esat ) + + + return end subroutine GetCanopyGasParameters + + ! ==================================================================================== + subroutine LeafLayerMaintenanceRespiration(lmr25top_ft, & + nscaler, & + ft, & + veg_tempk, & + lmr) + + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use pftconMod , only : pftcon + + ! Arguments + real(r8), intent(in) :: lmr25top_ft ! canopy top leaf maint resp rate at 25C + ! for this pft (umol CO2/m**2/s) + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) + + ! Locals + real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) + + + ! Parameter + real(r8), parameter :: lmrha = 46390._r8 ! activation energy for lmr (J/mol) + real(r8), parameter :: lmrhd = 150650._r8 ! deactivation energy for lmr (J/mol) + real(r8), parameter :: lmrse = 490._r8 ! entropy term for lmr (J/mol/K) + real(r8), parameter :: lmrc = 1.15912391_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + + ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s + ! ---------------------------------------------------------------------------------- + lmr25 = lmr25top_ft * nscaler + + if ( nint(pftcon%c3psn(ft)) == 1)then + lmr = lmr25 * ft1_f(veg_tempk, lmrha) * & + fth_f(veg_tempk, lmrhd, lmrse, lmrc) + else + lmr = lmr25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + lmr = lmr / (1._r8 + exp( 1.3_r8*(veg_tempk-(tfrz+55._r8)) )) + end if + + ! Any hydrodynamic limitations could go here, currently none + ! lmr = lmr * (nothing) + + end subroutine LeafLayerMaintenanceRespiration + + ! ==================================================================================== + + subroutine LeafLayerBiophysicalRates( parsun_lsl, & + ft, & + vcmax25top_ft, & + jmax25top_ft, & + tpu25top_ft, & + co2_rcurve_islope25top_ft, & + nscaler, & + veg_tempk, & + btran, & + vcmax, & + jmax, & + tpu, & + co2_rcurve_islope ) + + ! --------------------------------------------------------------------------------- + ! This subroutine calculates the localized rates of several key photosynthesis + ! rates. By localized, we mean specific to the plant type and leaf layer, + ! which factors in leaf physiology, as well as environmental effects. + ! This procedure should be called prior to iterative solvers, and should + ! have pre-calculated the reference rates for the pfts before this. + ! + ! The output biophysical rates are: + ! vcmax: maximum rate of carboxilation, + ! jmax: maximum electron transport rate, + ! tpu: triose phosphate utilization rate and + ! co2_rcurve_islope: initial slope of CO2 response curve (C4 plants) + ! --------------------------------------------------------------------------------- + + use pftconMod , only : pftcon + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + + ! Arguments + ! ------------------------------------------------------------------------------ + + real(r8), intent(in) :: parsun_lsl ! PAR absorbed in sunlit leaves for this layer + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: vcmax25top_ft ! canopy top maximum rate of carboxylation at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: jmax25top_ft ! canopy top maximum electron transport rate at 25C + ! for this pft (umol electrons/m**2/s) + real(r8), intent(in) :: tpu25top_ft ! canopy top triose phosphate utilization rate at 25C + ! for this pft (umol CO2/m**2/s) + real(r8), intent(in) :: co2_rcurve_islope25top_ft ! initial slope of CO2 response curve + ! (C4 plants) at 25C, canopy top, this pft + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) + + real(r8), intent(out) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(out) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(out) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) + real(r8), intent(out) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) + + ! Locals + ! ------------------------------------------------------------------------------- + real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C (umol CO2/m**2/s) + real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C (umol electrons/m**2/s) + real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) + real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C + + + ! Parameters + ! --------------------------------------------------------------------------------- + real(r8), parameter :: vcmaxha = 65330._r8 ! activation energy for vcmax (J/mol) + real(r8), parameter :: jmaxha = 43540._r8 ! activation energy for jmax (J/mol) + real(r8), parameter :: tpuha = 53100._r8 ! activation energy for tpu (J/mol) + real(r8), parameter :: vcmaxhd = 149250._r8 ! deactivation energy for vcmax (J/mol) + real(r8), parameter :: jmaxhd = 152040._r8 ! deactivation energy for jmax (J/mol) + real(r8), parameter :: tpuhd = 150650._r8 ! deactivation energy for tpu (J/mol) + real(r8), parameter :: vcmaxse = 485._r8 ! entropy term for vcmax (J/mol/K) + real(r8), parameter :: jmaxse = 495._r8 ! entropy term for jmax (J/mol/K) + real(r8), parameter :: tpuse = 490._r8 ! entropy term for tpu (J/mol/K) + real(r8), parameter :: vcmaxc = 1.1534040_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8), parameter :: jmaxc = 1.1657242_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8), parameter :: tpuc = 1.1591239_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + + if ( parsun_lsl <= 0._r8) then ! night time + vcmax = 0._r8 + jmax = 0._r8 + tpu = 0._r8 + co2_rcurve_islope = 0._r8 + else ! day time + vcmax25 = vcmax25top_ft * nscaler + jmax25 = jmax25top_ft * nscaler + tpu25 = tpu25top_ft * nscaler + co2_rcurve_islope25 = co2_rcurve_islope25top_ft * nscaler + + ! Adjust for temperature + vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc) + jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc) + tpu = tpu25 * ft1_f(veg_tempk, tpuha) * fth_f(veg_tempk, tpuhd, tpuse, tpuc) + + if (nint(pftcon%c3psn(ft)) /= 1) then + vcmax = vcmax25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) + vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-veg_tempk ) )) + vcmax = vcmax / (1._r8 + exp( 0.3_r8*(veg_tempk-(tfrz+40._r8)) )) + end if + co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) !q10 response of product limited psn. + end if + + ! Adjust for water limitations + vcmax = vcmax * btran + + return + end subroutine LeafLayerBiophysicalRates + + + end module EDPhotosynthesisMod From dae1b6a8b4fa2a98e30f752c9010fba32b427adc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 2 Dec 2016 00:05:39 -0800 Subject: [PATCH 03/15] Turns out sapwood in photosynthesis is sometimes out of phase with that which calculated in the dynamics time-step. I removed the check during photosynthesis. --- biogeochem/EDCohortDynamicsMod.F90 | 1 - biogeophys/EDPhotosynthesisMod.F90 | 33 +++++++++++++----------------- 2 files changed, 14 insertions(+), 20 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index cdca9ec65b..229f4ef95e 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -236,7 +236,6 @@ subroutine allocate_live_biomass(cc_p,mode) currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & currentcohort%laimemory)*leaf_frac - else ! Leaves are on (leaves_off_switch==1) !the purpose of this section is to figure out the root and stem biomass when the leaves are off diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index f1dfe07ae6..c0414b7897 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -238,7 +238,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! currentPatch%present(:,:) call UpdateCanopyNCanNRadPresent(currentPatch) - + ! Part IV. Identify some environmentally derived parameters: ! These quantities are biologically irrelevant ! Michaelis-Menten constant for CO2 (Pa) @@ -267,9 +267,10 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) do FT = 1,numpft_ed !calculate patch and pft specific properties at canopy top. + ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) lnc = 1._r8 / (slatop(FT) * leafcn(FT)) - + !at the moment in ED we assume that there is no active N cycle. This should change, of course. FIX(RF,032414) Sep2011. vcmax25top(FT) = fnitr(FT) !fudge - shortcut using fnitr as a proxy for vcmax... @@ -342,7 +343,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) btran_eff = currentPatch%btran_ft(ft) !! end if - vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) !vegetation area index. + ! Vegetation area index + vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) if (iv == 1) then laican = laican + 0.5_r8 * vai else @@ -352,14 +354,13 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Scale for leaf nitrogen profile nscaler = exp(-kn(FT) * laican) - + call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in nscaler, & ! in ft, & ! in bc_in(s)%t_veg_pa(ifp), & ! in lmr_z(CL,FT,iv)) ! out - call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(CL,FT,iv), & ! in ft, & ! in vcmax25top(ft), & ! in @@ -410,7 +411,6 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) enddo !PFT enddo !CL - !==============================================================================! ! Unpack fluxes from arrays into cohorts !==============================================================================! @@ -418,11 +418,11 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) call currentPatch%set_root_fraction(bc_in(s)%depth_gl) if(currentPatch%countcohorts > 0.0)then !avoid errors caused by empty patches - + currentCohort => currentPatch%tallest ! Cohort loop - + do while (associated(currentCohort)) ! Cohort loop - + if(currentCohort%n > 0._r8)then ! Zero cohort flux accumulators. @@ -435,6 +435,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Select canopy layer and PFT. FT = currentCohort%pft !are we going to have ftindex? CL = currentCohort%canopy_layer + !------------------------------------------------------------------------------ ! Accumulate fluxes over the sub-canopy layers of each cohort. !------------------------------------------------------------------------------ @@ -498,20 +499,14 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) leaf_frac = 1.0_r8/(currentCohort%canopy_trim + EDecophyscon%sapwood_ratio(currentCohort%pft) * & currentCohort%hite + pftcon%froot_leaf(currentCohort%pft)) - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! THIS CALCULATION SHOULD BE MOVED TO THE ALLOMETRY MODULE (RGK 10-8-2016) ! ------ IT ALSO SHOULD ALREADY HAVE BEEN CALCULATED RIGHT? ! ------ CHANGING TO A CHECK ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - if ( abs(currentCohort%bsw - (EDecophyscon%sapwood_ratio(currentCohort%pft) * & - currentCohort%hite * (currentCohort%balive + currentCohort%laimemory)*leaf_frac) ) & - > 1e-9 ) then - write(fates_log(),*) 'Sapwood biomass calculated during photosynthesis' - write(fates_log(),*) 'does not match what is contained in cohort%bsw' - write(fates_log(),*) 'which is the prognostic variable. Stopping.' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + 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 @@ -641,7 +636,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end if currentPatch => currentPatch%younger - + end do end do !site loop From e9b5fd21905a0740abb10a080308c989488e2be5 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 2 Dec 2016 12:36:42 -0800 Subject: [PATCH 04/15] Added some documentation to the photosynthesis refactoring. --- biogeophys/EDPhotosynthesisMod.F90 | 58 +++++++++++++++++------------- 1 file changed, 34 insertions(+), 24 deletions(-) diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index c0414b7897..66260d66d4 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -1,11 +1,10 @@ -module EDPhotosynthesisMod +module FATESPhotosynthesisMod - !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------------- ! !DESCRIPTION: - ! Calculates the photosynthetic fluxes for the ED model - ! This code is equivalent to the 'photosynthesis' subroutine in PhotosynthesisMod.F90. - ! We have split this out to reduce merge conflicts until we can pull out - ! common code used in both the ED and CLM versions. + ! Calculates the photosynthetic fluxes for the FATES model + ! This code is similar to and was originally based off of the 'photosynthesis' + ! subroutine in the CLM model. ! ! Parameter for activation and deactivation energies were taken from: ! Activation energy, from: @@ -18,7 +17,6 @@ module EDPhotosynthesisMod ! ------------------------------------------------------------------------------------ ! !USES: - ! use abortutils, only : endrun use FatesGlobals, only : fates_log @@ -309,17 +307,24 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end do !FT + ! ------------------------------------------------------------------------ + ! Part VI: Loop over all leaf layers. + ! The concept of leaf layers is a result of the radiative transfer scheme. + ! A leaf layer has uniform radiation environment. Leaf layers are a group + ! of vegetation surfaces (stems and leaves) which inhabit the same + ! canopy-layer "CL", have the same functional type "ft" and within those + ! two partitions are further partitioned into vertical layers where + ! downwelling radiation attenuates in order. + ! In this phase we loop over the leaf layers and calculate the + ! photosynthesis and respiration of the layer (since all biophysical + ! properties are homogeneous). After this step, we can loop through + ! our cohort list, associate each cohort with its list of leaf-layers + ! and transfer these quantities to the cohort. + ! With plant hydraulics, we must realize that photosynthesis and + ! respiration will be different for leaves of each cohort in the leaf + ! layers, as they will have there own hydraulic limitations. + ! ------------------------------------------------------------------------ - ! If we are using plant hydro-dynamics, then several photosynthesis - ! variables will be available at the cohort scale, and not the - ! pft scale. So here we split and use different looping structures - ! ------------------------------------------------------------------ - ! if ( use_fates_plant_hydro ) - - - !==============================================================================! - ! Calculate Nitrogen scaling factors and photosynthetic parameters. - !==============================================================================! do CL = 1, NCL_p do FT = 1,numpft_ed @@ -412,7 +417,10 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) enddo !CL !==============================================================================! - ! Unpack fluxes from arrays into cohorts + ! Part VII: Transfer leaf biophysical rates (like maintenance respiration, + ! carbon assimilation and conductance) that are defined by the leaf layer + ! (which is area independent, ie /m2) onto each cohort (where the rates become + ! per cohort, ie /individual). !==============================================================================! call currentPatch%set_root_fraction(bc_in(s)%depth_gl) @@ -499,11 +507,13 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) leaf_frac = 1.0_r8/(currentCohort%canopy_trim + EDecophyscon%sapwood_ratio(currentCohort%pft) * & currentCohort%hite + pftcon%froot_leaf(currentCohort%pft)) - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! THIS CALCULATION SHOULD BE MOVED TO THE ALLOMETRY MODULE (RGK 10-8-2016) - ! ------ IT ALSO SHOULD ALREADY HAVE BEEN CALCULATED RIGHT? - ! ------ CHANGING TO A CHECK - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! ------------------------------------------------------------------ + ! Part VIII: Calculate maintenance respiration in the sapwood and + ! fine root pools. + ! ------------------------------------------------------------------ + + currentCohort%bsw = EDecophyscon%sapwood_ratio(currentCohort%pft) * & currentCohort%hite * (currentCohort%balive + currentCohort%laimemory)*leaf_frac @@ -1471,4 +1481,4 @@ end subroutine LeafLayerBiophysicalRates -end module EDPhotosynthesisMod +end module FATESPhotosynthesisMod From 242f39c8c45f06f6e1724a527d357b9977e26d4d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 2 Dec 2016 17:48:39 -0800 Subject: [PATCH 05/15] Modified the looping structure in photosynthesis to accomodate leaf layer variability at the cohort and pft level. Also modularized the scaling of leaf-level fluxes back to cohorts. --- biogeophys/EDPhotosynthesisMod.F90 | 529 ++++++++++++++++------------- main/EDTypesMod.F90 | 10 + 2 files changed, 311 insertions(+), 228 deletions(-) diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index 66260d66d4..33d1a21066 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -22,13 +22,14 @@ module FATESPhotosynthesisMod use FatesGlobals, only : fates_log use FatesConstantsMod, only : r8 => fates_r8 use shr_log_mod , only : errMsg => shr_log_errMsg + use EDTypesMod, only : use_fates_plant_hydro implicit none private ! PUBLIC MEMBER FUNCTIONS: - public :: Photosynthesis_ED !ED specific photosynthesis routine + public :: FATESPhotosynthesis !ED specific photosynthesis routine character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -42,7 +43,7 @@ module FATESPhotosynthesisMod contains !--------------------------------------------------------- - subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) + subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! @@ -62,7 +63,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) use EDParamsMod , only : ED_val_ag_biomass use EDSharedParamsMod , only : EDParamsShareInst 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 @@ -100,36 +101,45 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! ----------------------------------------------------------------------------------- type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type), pointer :: currentCohort - integer , parameter :: psn_type = 2 !c3 or c4. + ! ----------------------------------------------------------------------------------- + ! These three arrays hold leaf-level biophysical rates that are calculated + ! in one loop and then sent to the cohorts in another loop. If hydraulics are + ! on, we calculate a unique solution for each level-cohort-layer combination. + ! If we are not using hydraulics, we calculate a unique solution for each + ! level-pft-layer combination. Thus the following three arrays are statically + ! allocated for the maximum space of the two cases (numCohortsPerPatch) + ! ----------------------------------------------------------------------------------- + + ! leaf maintenance (dark) respiration + real(r8) :: lmr_z(cp_nclmax,mxpft,cp_nlevcan) - real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) - ! either from cohort or patch-pft - ! - ! Leaf photosynthesis parameters - ! Note: None of these variables need to be an array. We put them - ! in arrays only to enable user debugging diagnostics - real(r8) :: vcmax_z(cp_nclmax,mxpft,cp_nlevcan) ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8) :: jmax_z(cp_nclmax,mxpft,cp_nlevcan) ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: tpu_z(cp_nclmax,mxpft,cp_nlevcan) ! triose phosphate utilization rate (umol CO2/m**2/s) - real(r8) :: kp_z(cp_nclmax,mxpft,cp_nlevcan) ! initial slope of CO2 response curve (C4 plants) - 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) :: anet_av(cp_nclmax,mxpft,cp_nlevcan) ! net leaf photosynthesis (umol CO2/m**2/s) - ! averaged over sun and shade leaves. + ! stomatal resistance s/m + real(r8) :: rs_z(cp_nclmax,mxpft,cp_nlevcan) + + ! net leaf photosynthesis averaged over sun and shade leaves. (umol CO2/m**2/s) + real(r8) :: anet_av(cp_nclmax,mxpft,cp_nlevcan) + + ! Mask used to determine which leaf-layer biophysical rates have been + ! used already + logical :: rate_mask(cp_nclmax,mxpft,cp_nlevcan) + + real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation (umol co2/m**2/s) + real(r8) :: jmax_z ! leaf layer maximum electron transport rate (umol electrons/m**2/s) + real(r8) :: tpu_z ! leaf layer triose phosphate utilization rate (umol CO2/m**2/s) + real(r8) :: kp_z ! leaf layer initial slope of CO2 response curve (C4 plants) real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) real(r8) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) ! --------------------------------------------------------------- ! TO-DO: bbbopt is slated to be transferred to the parameter file ! ---------------------------------------------------------------- real(r8) :: bbbopt(psn_type) ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s) real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) - real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation at 25C (umol CO2/m**2/s) real(r8) :: jmax25top(mxpft) ! canopy top: maximum electron transport rate at 25C (umol electrons/m**2/s) @@ -139,27 +149,22 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! Other integer :: cl,s,iv,j,ps,ft,ifp ! indices + integer :: nv ! number of leaf layers integer :: NCL_p ! number of canopy layers in patch real(r8) :: cf ! s m**2/umol -> s/m real(r8) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) real(r8) :: ceair ! vapor pressure of air, constrained (Pa) real(r8) :: nscaler ! leaf nitrogen scaling coefficient real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass -! 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. - real(r8) :: laifrac real(r8) :: tcsoi ! Temperature response function for root respiration. real(r8) :: tcwood ! Temperature response function for wood - -! real(r8) :: coarse_wood_frac ! amount of woody biomass that is coarse... real(r8) :: tree_area real(r8) :: gs_cohort real(r8) :: rscanopy real(r8) :: elai - 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) @@ -178,6 +183,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) real(r8),parameter :: base_mr_20 = 2.525e-6_r8 + + + associate( & c3psn => pftcon%c3psn , & slatop => pftcon%slatop , & ! specific leaf area at top of canopy, projected area basis [m^2/gC] @@ -307,6 +315,9 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) end do !FT + call currentPatch%set_root_fraction(bc_in(s)%depth_gl) + + ! ------------------------------------------------------------------------ ! Part VI: Loop over all leaf layers. ! The concept of leaf layers is a result of the radiative transfer scheme. @@ -324,49 +335,74 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) ! respiration will be different for leaves of each cohort in the leaf ! layers, as they will have there own hydraulic limitations. ! ------------------------------------------------------------------------ + rate_mask(:,:,:) = .false. - do CL = 1, NCL_p - do FT = 1,numpft_ed + if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches + + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) ! Cohort loop + + ! Identify the canopy layer (CL), functional type (FT) + ! and the leaf layer (IV) for this cohort + FT = currentCohort%pft + CL = currentCohort%canopy_layer if(currentPatch%present(CL,FT) == 1)then ! are there any leaves of this pft in this layer? - + if(CL==NCL_p)then !are we in the top canopy layer or a shaded layer? laican = 0._r8 else laican = sum(currentPatch%canopy_layer_lai(CL+1:NCL_p)) end if - - ! Loop through canopy layers (above snow). Respiration needs to be - ! calculated every timestep. Others are calculated only if daytime - do iv = 1, currentPatch%nrad(CL,FT) - - !! if (use_fates_plant_hydro) then - !! !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) - !! !! btran = currentCohort%btran(iv) - !! else - bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(FT), 1._r8) - btran_eff = currentPatch%btran_ft(ft) - !! end if - - ! Vegetation area index - vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) - if (iv == 1) then - laican = laican + 0.5_r8 * vai - else - laican = laican + 0.5_r8 * (currentPatch%elai_profile(CL,FT,iv-1)+ & - currentPatch%esai_profile(CL,FT,iv-1))+vai - end if + + ! Loop over sublayers, only calculate the leaf-layers biophysical rates + ! if this unique set has not been calculated + ! In non-hydraulic runs, many cohorts of the same pft may share the + ! same leaf layer and the properties will be the same. + ! We will ignore these + + do iv = 1,currentCohort%nv - ! Scale for leaf nitrogen profile - nscaler = exp(-kn(FT) * laican) - - call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in + ! ----------------------------------------------------------- + ! If we are doing plant hydro-dynamics (or any run-type + ! where cohorts may generate different photosynthetic rates + ! of other cohorts in the same canopy-pft-layer combo), + ! we re-calculate the leaf biophysical rates for the + ! cohort-layer combo of interest. + ! but in the vanilla case, we only re-calculate if it has + ! not been done yet. + ! ----------------------------------------------------------- + + if ( .not.rate_mask(ft,cl,iv) .or. use_fates_plant_hydro ) then + + if (use_fates_plant_hydro) then + !! !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) + !! !! btran = currentCohort%btran(iv) + else + bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(FT), 1._r8) + btran_eff = currentPatch%btran_ft(ft) + end if + + ! Vegetation area index + vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) + if (iv == 1) then + laican = laican + 0.5_r8 * vai + else + laican = laican + 0.5_r8 * (currentPatch%elai_profile(CL,FT,iv-1)+ & + currentPatch%esai_profile(CL,FT,iv-1))+vai + end if + + ! Scale for leaf nitrogen profile + nscaler = exp(-kn(FT) * laican) + + + call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in nscaler, & ! in ft, & ! in bc_in(s)%t_veg_pa(ifp), & ! in lmr_z(CL,FT,iv)) ! out - call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(CL,FT,iv), & ! in + call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(CL,FT,iv), & ! in ft, & ! in vcmax25top(ft), & ! in jmax25top(ft), & ! in @@ -375,12 +411,12 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) nscaler, & ! in bc_in(s)%t_veg_pa(ifp), & ! in btran_eff, & ! in - vcmax_z(CL,FT,iv), & ! out - jmax_z(CL,FT,iv), & ! out - tpu_z(CL,FT,iv), & ! out - kp_z(CL,FT,iv) ) ! out + vcmax_z, & ! out + jmax_z, & ! out + tpu_z, & ! out + kp_z ) ! out - call LeafLayerPhotosynthesis(currentPatch%f_sun(CL,FT,iv), & ! in + call LeafLayerPhotosynthesis(currentPatch%f_sun(CL,FT,iv), & ! in currentPatch%ed_parsun_z(CL,FT,iv), & ! in currentPatch%ed_parsha_z(CL,FT,iv), & ! in currentPatch%ed_laisun_z(CL,FT,iv), & ! in @@ -388,10 +424,10 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentPatch%canopy_area_profile(CL,FT,iv), & ! in ft, & ! in nscaler, & ! in - vcmax_z(CL,FT,iv), & ! in - jmax_z(CL,FT,iv), & ! in - tpu_z(CL,FT,iv), & ! in - kp_z(CL,FT,iv), & ! in + vcmax_z, & ! in + jmax_z, & ! in + tpu_z, & ! in + kp_z, & ! in bc_in(s)%t_veg_pa(ifp), & ! in bc_in(s)%esat_tv_pa(ifp), & ! in bc_in(s)%forc_pbot, & ! in @@ -409,115 +445,71 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentPatch%psn_z(cl,ft,iv), & ! out rs_z(CL,FT,iv), & ! out anet_av(CL,FT,iv)) ! out - - - end do ! iv - end if !present - enddo !PFT - enddo !CL - !==============================================================================! - ! Part VII: Transfer leaf biophysical rates (like maintenance respiration, - ! carbon assimilation and conductance) that are defined by the leaf layer - ! (which is area independent, ie /m2) onto each cohort (where the rates become - ! per cohort, ie /individual). - !==============================================================================! - - call currentPatch%set_root_fraction(bc_in(s)%depth_gl) - - if(currentPatch%countcohorts > 0.0)then !avoid errors caused by empty patches - - currentCohort => currentPatch%tallest ! Cohort loop - - do while (associated(currentCohort)) ! Cohort loop - - if(currentCohort%n > 0._r8)then - - ! Zero cohort flux accumulators. - currentCohort%npp_tstep = 0.0_r8 - currentCohort%resp_tstep = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%resp_m = 0.0_r8 - - ! Select canopy layer and PFT. - FT = currentCohort%pft !are we going to have ftindex? - CL = currentCohort%canopy_layer - - !------------------------------------------------------------------------------ - ! Accumulate fluxes over the sub-canopy layers of each cohort. - !------------------------------------------------------------------------------ - ! 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(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? - - 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%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) = anet_av(cl,ft,1:currentCohort%nv) * umolC_to_kgC * dtime + rate_mask(ft,cl,iv) = .true. + end if + end do + + ! Zero cohort flux accumulators. + currentCohort%npp_tstep = 0.0_r8 + currentCohort%resp_tstep = 0.0_r8 + currentCohort%gpp_tstep = 0.0_r8 + currentCohort%rdark = 0.0_r8 + currentCohort%resp_m = 0.0_r8 + currentCohort%ts_net_uptake = 0.0_r8 + + ! --------------------------------------------------------------- + ! Part VII: Transfer leaf flux rates (like maintenance respiration, + ! carbon assimilation and conductance) that are defined by the leaf layer + ! (which is area independent, ie /m2) onto each cohort (where the rates become + ! per cohort, ie /individual). Most likely a sum over layers. + ! --------------------------------------------------------------- + nv = currentCohort%nv + call ScaleLeafLayerFluxToCohort(nv, & !in + currentPatch%psn_z(cl,ft,1:nv), & !in + lmr_z(cl,ft,1:nv), & !in + rs_z(cl,ft,1:nv), & !in + anet_av(cl,ft,1:nv), & !in + currentPatch%elai_profile(cl,ft,1:nv), & !in + currentCohort%c_area, & !in + currentCohort%n, & !in + currentCohort%treelai, & !in + currentCohort%treesai, & !in + bc_in(s)%rb_pa(ifp), & !in + currentCohort%gscan, & ! out + currentCohort%gpp_tstep, & ! out + currentCohort%rdark) ! out + + ! Net Uptake does not need to be scaled, just transfer directly + currentCohort%ts_net_uptake(1:nv) = anet_av(cl,ft,1:nv) * umolC_to_kgC else - + + ! In this case, the cohort had no leaves, + ! so no productivity,conductance, transpiration uptake + ! or dark respiration + currentCohort%gpp_tstep = 0.0_r8 currentCohort%rdark = 0.0_r8 currentCohort%gscan = 0.0_r8 currentCohort%ts_net_uptake(:) = 0.0_r8 - - end if - - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 832 ', currentCohort%gpp_tstep - - laifrac = (currentCohort%treelai+currentCohort%treesai)-(currentCohort%nv-1)*dinc_ed - - gs_cohort = 1.0_r8/(rs_z(cl,ft,currentCohort%nv)+bc_in(s)%rb_pa(ifp))*laifrac*tree_area - currentCohort%gscan = currentCohort%gscan+gs_cohort - - if ( DEBUG ) then - 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(fates_log(),*) 'EDPhoto 843 ', currentCohort%rdark - - - currentCohort%rdark = currentCohort%rdark + lmr_z(cl,ft,currentCohort%nv) * & - currentPatch%elai_profile(cl,ft,currentCohort%nv) * laifrac * tree_area - - ! 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)) - + + end if ! if(currentPatch%present(CL,FT) == 1)then + ! ------------------------------------------------------------------ ! Part VIII: Calculate maintenance respiration in the sapwood and ! fine root pools. ! ------------------------------------------------------------------ - - + + 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%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 ! We are using the fine-root C:N ratio as an approximation for @@ -529,13 +521,13 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) live_croot_n = (1.0_r8-ED_val_ag_biomass) * currentCohort%bsw / & frootcn(currentCohort%pft) froot_n = currentCohort%br / frootcn(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. !------------------------------------------------------------------------------ - + ! Live stem MR (kgC/plant/s) (above ground sapwood) ! ------------------------------------------------------------------ if (woody(ft) == 1) then @@ -545,8 +537,8 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) else currentCohort%livestem_mr = 0._r8 end if - - + + ! Fine Root MR (kgC/plant/s) ! ------------------------------------------------------------------ currentCohort%froot_mr = 0._r8 @@ -555,7 +547,7 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) currentCohort%froot_mr = currentCohort%froot_mr + & froot_n * base_mr_20 * tcsoi * currentPatch%rootfr_ft(ft,j) enddo - + ! Coarse Root MR (kgC/plant/s) (below ground sapwood) ! ------------------------------------------------------------------ if (woody(ft) == 1) then @@ -570,90 +562,73 @@ subroutine Photosynthesis_ED (nsites, sites,bc_in,bc_out,dtime) 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(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 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%rdark - + ! convert from kgC/indiv/s to kgC/indiv/timestep - currentCohort%resp_m = currentCohort%resp_m * dtime - currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime + currentCohort%resp_m = currentCohort%resp_m * dtime + currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime + currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime + + + 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 currentCohort%npp_tstep = currentCohort%gpp_tstep - currentCohort%resp_tstep ! kgC/indiv/ts - - !------------------------------------------------------------------------------ - ! Remove whole plant respiration from net uptake. (kgC/indiv/ts) - if(currentCohort%treelai > 0._r8)then - ! do iv =1,currentCohort%NV - ! currentCohort%year_net_uptake(iv) = currentCohort%year_net_uptake(iv) - & - ! (timestep_secs*(currentCohort%livestem_mr + currentCohort%livecroot_mr & - ! minus contribution to whole plant respn. - ! + currentCohort%froot_mr))/(currentCohort%treelai*currentCohort%c_area/currentCohort%n) - ! enddo - else !lai<0 - currentCohort%gpp_tstep = 0._r8 - currentCohort%resp_m = 0._r8 - currentCohort%gscan = 0._r8 - end if - else !pft<0 n<0 - 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 - currentCohort%ts_net_uptake(1:currentCohort%nv) = 0._r8 - end if !pft<0 n<0 - - bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + currentCohort%gpp_tstep - bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + currentCohort%resp_m - ! accumulate cohort level canopy conductances over whole area before dividing by total area. - bc_out(s)%gccanopy_pa(ifp) = bc_out(s)%gccanopy_pa(ifp) + currentCohort%gscan * & - currentCohort%n /currentPatch%total_canopy_area - - currentCohort => currentCohort%shorter - - enddo ! end cohort loop. - end if !count_cohorts is more than zero. + + bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + currentCohort%gpp_tstep + bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + currentCohort%resp_m + ! accumulate cohort level canopy conductances over whole area before dividing by total area. + bc_out(s)%gccanopy_pa(ifp) = bc_out(s)%gccanopy_pa(ifp) + currentCohort%gscan * & + currentCohort%n /currentPatch%total_canopy_area + + currentCohort => currentCohort%shorter + + enddo ! end cohort loop. + end if !count_cohorts is more than zero. + + + elai = calc_areaindex(currentPatch,'elai') + + bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) / currentPatch%area + bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) / currentPatch%area + if(bc_out(s)%gccanopy_pa(ifp) > 1._r8/rsmax0 .and. elai > 0.0_r8)then + rscanopy = (1.0_r8/bc_out(s)%gccanopy_pa(ifp))-bc_in(s)%rb_pa(ifp)/elai ! this needs to be resistance per unit leaf area. + else + rscanopy = rsmax0 + 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/umol_per_mmol !convert into umol m-2 s-1 then mmol m-2 s-1. + end if + currentPatch => currentPatch%younger - elai = calc_areaindex(currentPatch,'elai') - - bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) / currentPatch%area - bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) / currentPatch%area - if(bc_out(s)%gccanopy_pa(ifp) > 1._r8/rsmax0 .and. elai > 0.0_r8)then - rscanopy = (1.0_r8/bc_out(s)%gccanopy_pa(ifp))-bc_in(s)%rb_pa(ifp)/elai ! this needs to be resistance per unit leaf area. - else - rscanopy = rsmax0 - 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/umol_per_mmol !convert into umol m-2 s-1 then mmol m-2 s-1. - end if - - currentPatch => currentPatch%younger + end do - end do + end do !site loop - end do !site loop - - end associate -end subroutine Photosynthesis_ED - + end associate + end subroutine FATESPhotosynthesis + ! ======================================================================================= subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in @@ -1018,7 +993,105 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in end associate return end subroutine LeafLayerPhotosynthesis - + +! ======================================================================================= + +subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv + psn_llz, & ! in %psn_z(cl,ft,1:currentCohort%nv) + lmr_llz, & ! in lmr_z(cl,ft,1:currentCohort%nv) + rs_llz, & ! in rs_z(cl,ft,1:currentCohort%nv) + anet_av_llz, & ! in anet_av_z(cl,ft,1:currentCohort%nv) + elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) + c_area, & ! in currentCohort%c_area + nplant, & ! in currentCohort%n + treelai, & ! in currentCohort%treelai + treesai, & ! in currentCohort%treesai + rb, & ! in bc_in(s)%rb_pa(ifp) + gscan, & ! out currentCohort%gscan + gpp, & ! out currentCohort%gpp_tstep + rdark) ! out currentCohort%rdark + + ! ------------------------------------------------------------------------------------ + ! This subroutine effectively integrates leaf carbon fluxes over the + ! leaf layers to give cohort totals. + ! Some arguments have the suffix "_llz". This indicates that the vector + ! is stratefied in the leaf-layer (ll) dimension, and is a portion of the calling + ! array which has the "_z" tag, thus "llz". + ! ------------------------------------------------------------------------------------ + + use FatesConstantsMod, only : umolC_to_kgC + use EDTypesMod, only : dinc_ed + + ! Arguments + integer, intent(in) :: nv ! number of active leaf layers + real(r8), intent(in) :: psn_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: lmr_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: rs_llz(nv) ! s/m + real(r8), intent(in) :: anet_av_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: elai_llz(nv) ! exposed LAI per layer + real(r8), intent(in) :: c_area ! crown area m2/m2 + real(r8), intent(in) :: nplant ! indiv/m2 + real(r8), intent(in) :: treelai ! m2/m2 + real(r8), intent(in) :: treesai ! m2/m2 + real(r8), intent(in) :: rb ! boundary layer resistance (s/m) + + real(r8), intent(out) :: gscan ! Canopy conductance of the cohort m/s + real(r8), intent(out) :: gpp ! GPP (kgC/indiv/s) + real(r8), intent(out) :: rdark ! Dark Leaf Respiration (kgC/indiv/s) + + ! GPP IN THIS SUBROUTINE IS A RATE. THE CALLING ARGUMENT IS GPP_TSTEP. AFTER THIS + ! CALL THE RATE WILL BE MULTIPLIED BY THE INTERVAL TO GIVE THE INTEGRATED QUANT. + + + ! Locals + real(r8) :: tree_area + real(r8) :: laifrac + + ! Convert from umolC/m2leaf/s to umolC/indiv/s ( x canopy area x 1m2 leaf area). + tree_area = c_area/nplant + + ! The routine is only called if there are leaves. If there are leaves, + ! there is at least 1 layer + + laifrac = (treelai+treesai)-dble(nv-1)*dinc_ed + + ! Canopy Conductance + gscan = 1.0_r8/(rs_llz(nv)+rb)*laifrac*tree_area + + ! GPP + gpp = psn_llz(nv) * elai_llz(nv) * laifrac * tree_area + + ! Dark respiration + rdark = lmr_llz(nv) * elai_llz(nv) * laifrac * tree_area + + ! If there is more than one layer, add the sum over the others + if ( nv>1 ) then + gpp = gpp + sum(psn_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area + rdark = rdark + sum(lmr_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area + gscan = gscan + sum((1.0_r8/(rs_llz(1:nv-1) + rb ))) * tree_area + end if + + ! Convert dark respiration and GPP from umol/plant/s to kgC/plant/s + ! Note that + + rdark = rdark * umolC_to_kgC + gpp = gpp * umolC_to_kgC + + + if ( DEBUG ) then + write(fates_log(),*) 'EDPhoto 816 ', gpp + write(fates_log(),*) 'EDPhoto 817 ', psn_llz(1:nv) + write(fates_log(),*) 'EDPhoto 820 ', nv + write(fates_log(),*) 'EDPhoto 821 ', elai_llz(1:nv) + write(fates_log(),*) 'EDPhoto 843 ', rdark + write(fates_log(),*) 'EDPhoto 871 ', laifrac + write(fates_log(),*) 'EDPhoto 872 ', tree_area + write(fates_log(),*) 'EDPhoto 873 ', nv + endif + + return +end subroutine ScaleLeafLayerFluxToCohort + ! ======================================================================================= function ft1_f(tl, ha) result(ans) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index d02891cb28..613bd9ef69 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -149,6 +149,16 @@ module EDTypesMod ! HLM will interpret that the value should not be included in the average. real(r8) :: cp_hio_ignore_val + + + ! Module switches (this will be read in one day) + ! This variable only exists now to serve as a place holder + !!!!!!!!!! THIS SHOULD NOT BE SET TO TRUE !!!!!!!!!!!!!!!!! + logical,parameter :: use_fates_plant_hydro = .false. + + + + !************************************ !** COHORT type structure ** !************************************ From 4073f576cd1e2e40e9f214072d63325f8334d6f2 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 4 Dec 2016 16:17:41 -0800 Subject: [PATCH 06/15] Photosynthesis refactors: mostly formatting changes in this commit group. --- ...esisMod.F90 => FatesPhotosynthesisMod.F90} | 815 +++++++++--------- 1 file changed, 426 insertions(+), 389 deletions(-) rename biogeophys/{EDPhotosynthesisMod.F90 => FatesPhotosynthesisMod.F90} (72%) diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/FatesPhotosynthesisMod.F90 similarity index 72% rename from biogeophys/EDPhotosynthesisMod.F90 rename to biogeophys/FatesPhotosynthesisMod.F90 index 33d1a21066..7a313a6581 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/FatesPhotosynthesisMod.F90 @@ -14,6 +14,8 @@ module FATESPhotosynthesisMod ! High temperature deactivation, from: ! Leuning (2002) Plant, Cell and Environment 25:1205-1210 ! The factor "c" scales the deactivation to a value of 1.0 at 25C + ! Photosynthesis and stomatal conductance parameters, from: + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 ! ------------------------------------------------------------------------------------ ! !USES: @@ -33,7 +35,7 @@ module FATESPhotosynthesisMod character(len=*), parameter, private :: sourcefile = & __FILE__ - !------------------------------------------------------------------------------ + !------------------------------------------------------------------------------------- ! maximum stomatal resistance [s/m] (used across several procedures) real(r8),parameter :: rsmax0 = 2.e4_r8 @@ -41,20 +43,21 @@ module FATESPhotosynthesisMod logical :: DEBUG = .false. contains - - !--------------------------------------------------------- - subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) - + + !-------------------------------------------------------------------------------------- + + subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) - ! + ! ----------------------------------------------------------------------------------- ! !DESCRIPTION: ! Leaf photosynthesis and stomatal conductance calculation as described by ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 and extended to ! a multi-layer canopy - ! + ! ----------------------------------------------------------------------------------- + + ! !USES: - use abortutils , only : endrun 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 @@ -62,33 +65,24 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) use EDParamsMod , only : ED_val_grperc use EDParamsMod , only : ED_val_ag_biomass use EDSharedParamsMod , only : EDParamsShareInst - use EDTypesMod , only : numpft_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 use FatesInterfaceMod , only : bc_out_type - use EDCanopyStructureMod, only : calc_areaindex - 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: + ! ARGUMENTS: ! ----------------------------------------------------------------------------------- integer,intent(in) :: nsites type(ed_site_type),intent(inout),target :: sites(nsites) @@ -97,11 +91,11 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) real(r8),intent(in) :: dtime - ! !LOCAL VARIABLES: + ! LOCAL VARIABLES: ! ----------------------------------------------------------------------------------- type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type), pointer :: currentCohort - integer , parameter :: psn_type = 2 !c3 or c4. + ! ----------------------------------------------------------------------------------- ! These three arrays hold leaf-level biophysical rates that are calculated ! in one loop and then sent to the cohorts in another loop. If hydraulics are @@ -111,18 +105,18 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! allocated for the maximum space of the two cases (numCohortsPerPatch) ! ----------------------------------------------------------------------------------- - ! leaf maintenance (dark) respiration + ! leaf maintenance (dark) respiration (umol CO2/m**2/s) Double check this real(r8) :: lmr_z(cp_nclmax,mxpft,cp_nlevcan) ! stomatal resistance s/m real(r8) :: rs_z(cp_nclmax,mxpft,cp_nlevcan) ! net leaf photosynthesis averaged over sun and shade leaves. (umol CO2/m**2/s) - real(r8) :: anet_av(cp_nclmax,mxpft,cp_nlevcan) + real(r8) :: anet_av_z(cp_nclmax,mxpft,cp_nlevcan) ! Mask used to determine which leaf-layer biophysical rates have been ! used already - logical :: rate_mask(cp_nclmax,mxpft,cp_nlevcan) + logical :: rate_mask_z(cp_nclmax,mxpft,cp_nlevcan) real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation (umol co2/m**2/s) @@ -135,10 +129,7 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) real(r8) :: co2_cpoint ! CO2 compensation point (Pa) real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) - ! --------------------------------------------------------------- - ! TO-DO: bbbopt is slated to be transferred to the parameter file - ! ---------------------------------------------------------------- - real(r8) :: bbbopt(psn_type) ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s) + real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation at 25C (umol CO2/m**2/s) @@ -183,7 +174,17 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) real(r8),parameter :: base_mr_20 = 2.525e-6_r8 - + ! ----------------------------------------------------------------------------------- + ! Photosynthesis and stomatal conductance parameters, from: + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 + ! ----------------------------------------------------------------------------------- + + ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s) + ! For C3 and C4 plants + ! ----------------------------------------------------------------------------------- + ! TO-DO: bbbopt is slated to be transferred to the parameter file + ! ----------------------------------------------------------------------------------- + real(r8),parameter, dimension(2) :: bbbopt = [10000._r8,40000._r8] associate( & @@ -196,17 +197,6 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) frootcn => pftcon%frootcn , & ! froot C:N (gc/gN) ! slope of BB relationship q10 => EDParamsShareInst%Q10) - - !==============================================================================! - ! Photosynthesis and stomatal conductance parameters, from: - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - !==============================================================================! - - ! Miscellaneous parameters, from Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - - bbbopt(1) = 10000._r8 - bbbopt(2) = 40000._r8 - do s = 1,nsites ! Multi-layer parameters scaled by leaf nitrogen profile. @@ -238,8 +228,9 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) if(bc_in(s)%filter_photo_pa(ifp)==2)then - ! Part III. Calculate the number of sublayers for each pft and layer. And then identify - ! which layer/pft combinations have things in them. Output: + ! Part III. Calculate the number of sublayers for each pft and layer. + ! And then identify which layer/pft combinations have things in them. + ! Output: ! currentPatch%ncan(:,:) ! currentPatch%present(:,:) call UpdateCanopyNCanNRadPresent(currentPatch) @@ -271,34 +262,42 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! but not environmentally dependent ! ------------------------------------------------------------------------ - do FT = 1,numpft_ed !calculate patch and pft specific properties at canopy top. + do ft = 1,numpft_ed ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) - lnc = 1._r8 / (slatop(FT) * leafcn(FT)) + lnc = 1._r8 / (slatop(ft) * leafcn(ft)) - !at the moment in ED we assume that there is no active N cycle. This should change, of course. FIX(RF,032414) Sep2011. - vcmax25top(FT) = fnitr(FT) !fudge - shortcut using fnitr as a proxy for vcmax... + ! at the moment in ED we assume that there is no active N cycle. + ! This should change, of course. FIX(RF,032414) Sep2011. + ! fudge - shortcut using fnitr as a proxy for vcmax... + vcmax25top(ft) = fnitr(ft) - ! Parameters derived from vcmax25top. Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of Experimental Botany 44:907-920. - ! Here use a factor "1.67", from Medlyn et al (2002) Plant, Cell and Environment 25:1167-1179 + ! Parameters derived from vcmax25top. + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 + ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of + ! Experimental Botany 44:907-920. 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)-tfrzc),11._r8),35._r8)) * vcmax25top(FT) + ! 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)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) - jmax25top(FT) = 1.67_r8 * vcmax25top(FT) - tpu25top(FT) = 0.167_r8 * vcmax25top(FT) - kp25top(FT) = 20000._r8 * vcmax25top(FT) - - ! Nitrogen scaling factor. Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used - ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al (2010) Biogeosciences, 7, 1833-1859 + jmax25top(ft) = 1.67_r8 * vcmax25top(ft) + tpu25top(ft) = 0.167_r8 * vcmax25top(ft) + kp25top(ft) = 20000._r8 * vcmax25top(ft) + + ! Nitrogen scaling factor. + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used + ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al + ! (2010) Biogeosciences, 7, 1833-1859 ! Remove daylength factor from vcmax25 so that kn is based on maximum vcmax25 if (bc_in(s)%dayl_factor_pa(ifp) == 0._r8) then - kn(FT) = 0._r8 + kn(ft) = 0._r8 else - kn(FT) = exp(0.00963_r8 * vcmax25top(FT) - 2.43_r8) + kn(ft) = exp(0.00963_r8 * vcmax25top(ft) - 2.43_r8) end if ! Leaf maintenance respiration to match the base rate used in CN @@ -310,14 +309,13 @@ subroutine FATESPhotosynthesis (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 / (umolC_to_kgC * g_per_kg) + lmr25top(ft) = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) + lmr25top(ft) = lmr25top(ft) * lnc / (umolC_to_kgC * g_per_kg) - end do !FT + end do !ft call currentPatch%set_root_fraction(bc_in(s)%depth_gl) - ! ------------------------------------------------------------------------ ! Part VI: Loop over all leaf layers. ! The concept of leaf layers is a result of the radiative transfer scheme. @@ -335,35 +333,32 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! respiration will be different for leaves of each cohort in the leaf ! layers, as they will have there own hydraulic limitations. ! ------------------------------------------------------------------------ - rate_mask(:,:,:) = .false. + rate_mask_z(:,:,:) = .false. if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches currentCohort => currentPatch%tallest do while (associated(currentCohort)) ! Cohort loop - ! Identify the canopy layer (CL), functional type (FT) + ! Identify the canopy layer (cl), functional type (ft) ! and the leaf layer (IV) for this cohort - FT = currentCohort%pft - CL = currentCohort%canopy_layer + ft = currentCohort%pft + cl = currentCohort%canopy_layer - if(currentPatch%present(CL,FT) == 1)then ! are there any leaves of this pft in this layer? + + ! are there any leaves of this pft in this layer? + if(currentPatch%present(cl,ft) == 1)then - if(CL==NCL_p)then !are we in the top canopy layer or a shaded layer? + if(cl==NCL_p)then !are we in the top canopy layer or a shaded layer? laican = 0._r8 else - laican = sum(currentPatch%canopy_layer_lai(CL+1:NCL_p)) + laican = sum(currentPatch%canopy_layer_lai(cl+1:NCL_p)) end if - - ! Loop over sublayers, only calculate the leaf-layers biophysical rates - ! if this unique set has not been calculated - ! In non-hydraulic runs, many cohorts of the same pft may share the - ! same leaf layer and the properties will be the same. - ! We will ignore these + ! Loop over leaf-layers do iv = 1,currentCohort%nv - ! ----------------------------------------------------------- + ! ------------------------------------------------------------ ! If we are doing plant hydro-dynamics (or any run-type ! where cohorts may generate different photosynthetic rates ! of other cohorts in the same canopy-pft-layer combo), @@ -371,82 +366,100 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! cohort-layer combo of interest. ! but in the vanilla case, we only re-calculate if it has ! not been done yet. - ! ----------------------------------------------------------- - - if ( .not.rate_mask(ft,cl,iv) .or. use_fates_plant_hydro ) then + ! ------------------------------------------------------------ + + if ( .not.rate_mask_z(ft,cl,iv) .or. use_fates_plant_hydro ) then if (use_fates_plant_hydro) then + write(fates_log(),*) 'use_fates_plant_hydro in EDTypes' + write(fates_log(),*) 'has been set to true. You have inadvartently' + write(fates_log(),*) 'turned on a future feature that is not in the' + write(fates_log(),*) 'FATES model codeset yet. Please set this to' + write(fates_log(),*) 'false and re-compile.' + call endrun(msg=errMsg(sourcefile, __LINE__)) !! !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) !! !! btran = currentCohort%btran(iv) else - bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(FT), 1._r8) + bbb = max (bbbopt(nint(c3psn(ft)))*currentPatch%btran_ft(ft), 1._r8) btran_eff = currentPatch%btran_ft(ft) end if - + ! Vegetation area index - vai = (currentPatch%elai_profile(CL,FT,iv)+currentPatch%esai_profile(CL,FT,iv)) + vai = (currentPatch%elai_profile(cl,ft,iv)+currentPatch%esai_profile(cl,ft,iv)) if (iv == 1) then laican = laican + 0.5_r8 * vai else - laican = laican + 0.5_r8 * (currentPatch%elai_profile(CL,FT,iv-1)+ & - currentPatch%esai_profile(CL,FT,iv-1))+vai + laican = laican + 0.5_r8 * (currentPatch%elai_profile(cl,ft,iv-1)+ & + currentPatch%esai_profile(cl,ft,iv-1))+vai end if ! Scale for leaf nitrogen profile - nscaler = exp(-kn(FT) * laican) - + nscaler = exp(-kn(ft) * laican) + ! Part VII: Calculate dark respiration (leaf maintenance) for this layer call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in - nscaler, & ! in - ft, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - lmr_z(CL,FT,iv)) ! out - - call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(CL,FT,iv), & ! in - ft, & ! in - vcmax25top(ft), & ! in - jmax25top(ft), & ! in - tpu25top(ft), & ! in - kp25top(ft), & ! in - nscaler, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - btran_eff, & ! in - vcmax_z, & ! out - jmax_z, & ! out - tpu_z, & ! out - kp_z ) ! out - - call LeafLayerPhotosynthesis(currentPatch%f_sun(CL,FT,iv), & ! in - currentPatch%ed_parsun_z(CL,FT,iv), & ! in - currentPatch%ed_parsha_z(CL,FT,iv), & ! in - currentPatch%ed_laisun_z(CL,FT,iv), & ! in - currentPatch%ed_laisha_z(CL,FT,iv), & ! in - currentPatch%canopy_area_profile(CL,FT,iv), & ! in - ft, & ! in - nscaler, & ! in - vcmax_z, & ! in - jmax_z, & ! in - tpu_z, & ! in - kp_z, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - bc_in(s)%forc_pbot, & ! in - bc_in(s)%cair_pa(ifp), & ! in - bc_in(s)%oair_pa(ifp), & ! in - btran_eff, & ! in - bbb, & ! in - cf, & ! in - gb_mol, & ! in - ceair, & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr_z(CL,FT,iv), & ! in - currentPatch%psn_z(cl,ft,iv), & ! out - rs_z(CL,FT,iv), & ! out - anet_av(CL,FT,iv)) ! out - - rate_mask(ft,cl,iv) = .true. + nscaler, & ! in + ft, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + lmr_z(cl,ft,iv)) ! out + + ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), + ! (2) maximum electron transport rate, (3) triose phosphate + ! utilization rate and (4) the initial slope of CO2 response curve + ! (C4 plants). Earlier we calculated their base rates as dictated + ! by their plant functional type and some simple scaling rules for + ! nitrogen limitation baesd on canopy position (not prognostic). + ! These rates are the specific rates used in the actual photosynthesis + ! calculations that take localized environmental effects (temperature) + ! into consideration. + + call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(cl,ft,iv), & ! in + ft, & ! in + vcmax25top(ft), & ! in + jmax25top(ft), & ! in + tpu25top(ft), & ! in + kp25top(ft), & ! in + nscaler, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + btran_eff, & ! in + vcmax_z, & ! out + jmax_z, & ! out + tpu_z, & ! out + kp_z ) ! out + + ! Part IX: This call calculates the actual photosynthesis for the + ! leaf layer, as well as the stomatal resistance and the net assimilated carbon. + + call LeafLayerPhotosynthesis(currentPatch%f_sun(cl,ft,iv), & ! in + currentPatch%ed_parsun_z(cl,ft,iv), & ! in + currentPatch%ed_parsha_z(cl,ft,iv), & ! in + currentPatch%ed_laisun_z(cl,ft,iv), & ! in + currentPatch%ed_laisha_z(cl,ft,iv), & ! in + currentPatch%canopy_area_profile(cl,ft,iv), & ! in + ft, & ! in + vcmax_z, & ! in + jmax_z, & ! in + tpu_z, & ! in + kp_z, & ! in + bc_in(s)%t_veg_pa(ifp), & ! in + bc_in(s)%esat_tv_pa(ifp), & ! in + bc_in(s)%forc_pbot, & ! in + bc_in(s)%cair_pa(ifp), & ! in + bc_in(s)%oair_pa(ifp), & ! in + btran_eff, & ! in + bbb, & ! in + cf, & ! in + gb_mol, & ! in + ceair, & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr_z(cl,ft,iv), & ! in + currentPatch%psn_z(cl,ft,iv), & ! out + rs_z(cl,ft,iv), & ! out + anet_av_z(cl,ft,iv)) ! out + + rate_mask_z(ft,cl,iv) = .true. end if end do @@ -460,28 +473,29 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! --------------------------------------------------------------- ! Part VII: Transfer leaf flux rates (like maintenance respiration, - ! carbon assimilation and conductance) that are defined by the leaf layer - ! (which is area independent, ie /m2) onto each cohort (where the rates become - ! per cohort, ie /individual). Most likely a sum over layers. + ! carbon assimilation and conductance) that are defined by the + ! leaf layer (which is area independent, ie /m2) onto each cohort + ! (where the rates become per cohort, ie /individual). Most likely + ! a sum over layers. ! --------------------------------------------------------------- nv = currentCohort%nv call ScaleLeafLayerFluxToCohort(nv, & !in currentPatch%psn_z(cl,ft,1:nv), & !in lmr_z(cl,ft,1:nv), & !in rs_z(cl,ft,1:nv), & !in - anet_av(cl,ft,1:nv), & !in + anet_av_z(cl,ft,1:nv), & !in currentPatch%elai_profile(cl,ft,1:nv), & !in currentCohort%c_area, & !in currentCohort%n, & !in currentCohort%treelai, & !in currentCohort%treesai, & !in bc_in(s)%rb_pa(ifp), & !in - currentCohort%gscan, & ! out - currentCohort%gpp_tstep, & ! out - currentCohort%rdark) ! out + currentCohort%gscan, & !out + currentCohort%gpp_tstep, & !out + currentCohort%rdark) !out ! Net Uptake does not need to be scaled, just transfer directly - currentCohort%ts_net_uptake(1:nv) = anet_av(cl,ft,1:nv) * umolC_to_kgC + currentCohort%ts_net_uptake(1:nv) = anet_av_z(cl,ft,1:nv) * umolC_to_kgC else @@ -494,20 +508,22 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) currentCohort%gscan = 0.0_r8 currentCohort%ts_net_uptake(:) = 0.0_r8 - end if ! if(currentPatch%present(CL,FT) == 1)then + end if ! if(currentPatch%present(cl,ft) == 1)then ! ------------------------------------------------------------------ ! Part VIII: Calculate maintenance respiration in the sapwood and ! fine root pools. ! ------------------------------------------------------------------ - - leaf_frac = 1.0_r8/(currentCohort%canopy_trim + EDecophyscon%sapwood_ratio(currentCohort%pft) * & - currentCohort%hite + pftcon%froot_leaf(currentCohort%pft)) + + 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%hite * & + (currentCohort%balive + currentCohort%laimemory)*leaf_frac ! Calculate the amount of nitrogen in the above and below ground @@ -524,8 +540,12 @@ subroutine FATESPhotosynthesis (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. + ! Calculate Whole Plant Respiration + ! (this doesn't really need to be in this iteration at all, surely?) + ! Response: (RGK 12-2016): I think the positioning of these calls is + ! appropriate as of now. Maintenance calculations in sapwood and roots + ! vary by cohort and with changing temperature at the minimum, and there are + ! no sub-pools chopping up those pools any finer that need to be dealty with. !------------------------------------------------------------------------------ ! Live stem MR (kgC/plant/s) (above ground sapwood) @@ -574,31 +594,43 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! 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%livestem_mr + & + currentCohort%livecroot_mr + & + currentCohort%froot_mr + + ! no drought response right now.. something like: + ! resp_m = resp_m * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft) * & + ! pftcon%resp_drought_response(ft)) + currentCohort%resp_m = currentCohort%resp_m + currentCohort%rdark ! convert from kgC/indiv/s to kgC/indiv/timestep - currentCohort%resp_m = currentCohort%resp_m * dtime - currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime - currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime - - - + currentCohort%resp_m = currentCohort%resp_m * dtime + currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime + currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime 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 - currentCohort%npp_tstep = currentCohort%gpp_tstep - currentCohort%resp_tstep ! kgC/indiv/ts + 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 - bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + currentCohort%gpp_tstep - bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + currentCohort%resp_m - ! accumulate cohort level canopy conductances over whole area before dividing by total area. - bc_out(s)%gccanopy_pa(ifp) = bc_out(s)%gccanopy_pa(ifp) + currentCohort%gscan * & - currentCohort%n /currentPatch%total_canopy_area + bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + & + currentCohort%gpp_tstep + bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + & + currentCohort%resp_m + + ! accumulate cohort level canopy conductances over + ! whole area before dividing by total area + bc_out(s)%gccanopy_pa(ifp) = bc_out(s)%gccanopy_pa(ifp) + & + currentCohort%gscan * & + currentCohort%n /currentPatch%total_canopy_area currentCohort => currentCohort%shorter @@ -610,14 +642,16 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) / currentPatch%area bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) / currentPatch%area + if(bc_out(s)%gccanopy_pa(ifp) > 1._r8/rsmax0 .and. elai > 0.0_r8)then - rscanopy = (1.0_r8/bc_out(s)%gccanopy_pa(ifp))-bc_in(s)%rb_pa(ifp)/elai ! this needs to be resistance per unit leaf area. + rscanopy = (1.0_r8/bc_out(s)%gccanopy_pa(ifp))-bc_in(s)%rb_pa(ifp)/elai else rscanopy = rsmax0 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/umol_per_mmol !convert into umol m-2 s-1 then mmol m-2 s-1. + !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 end if currentPatch => currentPatch%younger @@ -627,78 +661,76 @@ subroutine FATESPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) end do !site loop end associate - end subroutine FATESPhotosynthesis - -! ======================================================================================= - -subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in - parsun_lsl, & ! in - parsha_lsl, & ! in - laisun_lsl, & ! in - laisha_lsl, & ! in - canopy_area_lsl, & ! in - ft, & ! in - nscaler, & ! in - vcmax, & ! in - jmax, & ! in - tpu, & ! in - co2_rcurve_islope, & ! in - veg_tempk, & ! in - veg_esat, & ! in - can_press, & ! in - can_co2_ppress, & ! in - can_o2_ppress, & ! in - btran, & ! in - bbb, & ! in - cf, & ! in - gb_mol, & ! in - ceair, & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr, & ! in - psn_out, & ! out - rstoma_out, & ! out - anet_av_out) ! out - - ! ------------------------------------------------------------------------------------ - ! This subroutine calculates photosynthesis and stomatal conductance within each leaf - ! sublayer. - ! A note on naming conventions: As this subroutine is called for every - ! leaf-sublayer, many of the arguments are specific to that "leaf sub layer" - ! (LSL), those variables are given a dimension tag "_lsl" - ! Other arguments or variables may be indicative of scales broader than the LSL. - ! ------------------------------------------------------------------------------------ - - use EDEcophysContype , only : EDecophyscon - use pftconMod , only : pftcon - - ! Arguments - ! ------------------------------------------------------------------------ - real(r8), intent(in) :: f_sun_lsl - real(r8), intent(in) :: parsun_lsl - real(r8), intent(in) :: parsha_lsl - real(r8), intent(in) :: laisun_lsl - real(r8), intent(in) :: laisha_lsl - real(r8), intent(in) :: canopy_area_lsl - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8), intent(in) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) - real(r8), intent(in) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) - real(r8), intent(in) :: veg_tempk ! vegetation temperature - real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) - - ! Important Note on the following gas pressures. This photosynthesis scheme will iteratively - ! solve for the co2 partial pressure at the leaf surface (ie in the stomata). The reference - ! point for these input values are NOT within that boundary layer that separates the stomata from - ! the canopy air space. The reference point for these is on the outside of that boundary - ! layer. This routine, which operates at the leaf scale, makes no assumptions about what the - ! scale of the refernce is, it could be lower atmosphere, it could be within the canopy - ! but most likely it is the closest value one can get to the edge of the leaf's boundary - ! layer. We use the convention "can_" because a reference point of within the canopy - ! ia a best reasonable scenario of where we can get that information from. + end subroutine FatesPhotosynthesis + + ! ======================================================================================= + + subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in + parsun_lsl, & ! in + parsha_lsl, & ! in + laisun_lsl, & ! in + laisha_lsl, & ! in + canopy_area_lsl, & ! in + ft, & ! in + vcmax, & ! in + jmax, & ! in + tpu, & ! in + co2_rcurve_islope, & ! in + veg_tempk, & ! in + veg_esat, & ! in + can_press, & ! in + can_co2_ppress, & ! in + can_o2_ppress, & ! in + btran, & ! in + bbb, & ! in + cf, & ! in + gb_mol, & ! in + ceair, & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr, & ! in + psn_out, & ! out + rstoma_out, & ! out + anet_av_out) ! out + + ! ------------------------------------------------------------------------------------ + ! This subroutine calculates photosynthesis and stomatal conductance within each leaf + ! sublayer. + ! A note on naming conventions: As this subroutine is called for every + ! leaf-sublayer, many of the arguments are specific to that "leaf sub layer" + ! (LSL), those variables are given a dimension tag "_lsl" + ! Other arguments or variables may be indicative of scales broader than the LSL. + ! ------------------------------------------------------------------------------------ + + use EDEcophysContype , only : EDecophyscon + use pftconMod , only : pftcon + + ! Arguments + ! ------------------------------------------------------------------------------------ + real(r8), intent(in) :: f_sun_lsl ! + real(r8), intent(in) :: parsun_lsl ! Absorbed PAR in sunlist leaves + real(r8), intent(in) :: parsha_lsl ! Absorved PAR in shaded leaves + real(r8), intent(in) :: laisun_lsl ! LAI in sunlit leaves + real(r8), intent(in) :: laisha_lsl ! LAI in shaded leaves + real(r8), intent(in) :: canopy_area_lsl ! + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8), intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8), intent(in) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) + real(r8), intent(in) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) + + ! Important Note on the following gas pressures. This photosynthesis scheme will iteratively + ! solve for the co2 partial pressure at the leaf surface (ie in the stomata). The reference + ! point for these input values are NOT within that boundary layer that separates the stomata from + ! the canopy air space. The reference point for these is on the outside of that boundary + ! layer. This routine, which operates at the leaf scale, makes no assumptions about what the + ! scale of the refernce is, it could be lower atmosphere, it could be within the canopy + ! but most likely it is the closest value one can get to the edge of the leaf's boundary + ! layer. We use the convention "can_" because a reference point of within the canopy + ! ia a best reasonable scenario of where we can get that information from. real(r8), intent(in) :: can_press ! Air pressure NEAR the surface of the leaf (Pa) real(r8), intent(in) :: can_co2_ppress ! Partial pressure of CO2 NEAR the leaf surface (Pa) @@ -721,7 +753,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Locals ! ------------------------------------------------------------------------ - integer :: ps ! Index for the different photosynthetic pathways C3,C4 + integer :: pp_type ! Index for the different photosynthetic pathways C3,C4 integer :: sunsha ! Index for differentiating sun and shade real(r8) :: gstoma ! Stomatal Conductance of this leaf layer (m/s) real(r8) :: agross ! co-limited gross leaf photosynthesis (umol CO2/m**2/s) @@ -765,14 +797,12 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! empirical curvature parameter for ap photosynthesis co-limitation real(r8),parameter :: theta_ip = 0.95_r8 - - associate( c3psn => pftcon%c3psn, & ! photosynthetic pathway: 0. = c4, 1. = c3 - bb_slope => EDecophyscon%BB_slope ) ! slope of BB relationship + associate( bb_slope => EDecophyscon%BB_slope ) ! slope of BB relationship - if (nint(c3psn(ft)) == 1)then - ps = 1 + if (nint(pftcon%c3psn(ft)) == 1) then! photosynthetic pathway: 0. = c4, 1. = c3 + pp_type = 1 else - ps = 2 + pp_type = 2 end if ! Part III: Photosynthesis and Conductance @@ -832,7 +862,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! Iterative loop for ci beginning with initial guess ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) - if (nint(c3psn(FT)) == 1)then + if (pp_type == 1)then co2_intra_c = init_a2l_co2_c3 * can_co2_ppress else co2_intra_c = init_a2l_co2_c4 * can_co2_ppress @@ -848,7 +878,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in co2_intra_c_old = co2_intra_c ! Photosynthesis limitation rate calculations - if (nint(c3psn(FT)) == 1)then + if (pp_type == 1)then ! C3: Rubisco-limited photosynthesis ac = vcmax * max(co2_intra_c-co2_cpoint, 0._r8) / & @@ -868,14 +898,14 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! C4: RuBP-limited photosynthesis if(sunsha == 1)then !sunlit if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then !guard against /0's in the night. - aj = quant_eff(ps) * parsun_lsl * 4.6_r8 + aj = quant_eff(pp_type) * parsun_lsl * 4.6_r8 !convert from per cohort to per m2 of leaf) aj = aj / (laisun_lsl * canopy_area_lsl) else aj = 0._r8 end if else - aj = quant_eff(ps) * parsha_lsl * 4.6_r8 + aj = quant_eff(pp_type) * parsha_lsl * 4.6_r8 aj = aj / (laisha_lsl * canopy_area_lsl) end if @@ -885,7 +915,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in end if ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap - aquad = theta_cj(ps) + aquad = theta_cj(pp_type) bquad = -(ac + aj) cquad = ac * aj call quadratic_f (aquad, bquad, cquad, r1, r2) @@ -976,9 +1006,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in write (fates_log(),*) 'CF: Ball-Berry error check - stomatal conductance error:' write (fates_log(),*) gs_mol, gs_mol_err end if - + enddo !sunsha loop - + !average leaf-level stomatal resistance rate over sun and shade leaves... rstoma_out = 1._r8/gstoma @@ -992,111 +1022,108 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in end if ! night or day end associate return -end subroutine LeafLayerPhotosynthesis - -! ======================================================================================= - -subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv - psn_llz, & ! in %psn_z(cl,ft,1:currentCohort%nv) - lmr_llz, & ! in lmr_z(cl,ft,1:currentCohort%nv) - rs_llz, & ! in rs_z(cl,ft,1:currentCohort%nv) - anet_av_llz, & ! in anet_av_z(cl,ft,1:currentCohort%nv) - elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) - c_area, & ! in currentCohort%c_area - nplant, & ! in currentCohort%n - treelai, & ! in currentCohort%treelai - treesai, & ! in currentCohort%treesai - rb, & ! in bc_in(s)%rb_pa(ifp) - gscan, & ! out currentCohort%gscan - gpp, & ! out currentCohort%gpp_tstep - rdark) ! out currentCohort%rdark - - ! ------------------------------------------------------------------------------------ - ! This subroutine effectively integrates leaf carbon fluxes over the - ! leaf layers to give cohort totals. - ! Some arguments have the suffix "_llz". This indicates that the vector - ! is stratefied in the leaf-layer (ll) dimension, and is a portion of the calling - ! array which has the "_z" tag, thus "llz". - ! ------------------------------------------------------------------------------------ - - use FatesConstantsMod, only : umolC_to_kgC - use EDTypesMod, only : dinc_ed - - ! Arguments - integer, intent(in) :: nv ! number of active leaf layers - real(r8), intent(in) :: psn_llz(nv) ! umolC/m2leaf/s - real(r8), intent(in) :: lmr_llz(nv) ! umolC/m2leaf/s - real(r8), intent(in) :: rs_llz(nv) ! s/m - real(r8), intent(in) :: anet_av_llz(nv) ! umolC/m2leaf/s - real(r8), intent(in) :: elai_llz(nv) ! exposed LAI per layer - real(r8), intent(in) :: c_area ! crown area m2/m2 - real(r8), intent(in) :: nplant ! indiv/m2 - real(r8), intent(in) :: treelai ! m2/m2 - real(r8), intent(in) :: treesai ! m2/m2 - real(r8), intent(in) :: rb ! boundary layer resistance (s/m) - - real(r8), intent(out) :: gscan ! Canopy conductance of the cohort m/s - real(r8), intent(out) :: gpp ! GPP (kgC/indiv/s) - real(r8), intent(out) :: rdark ! Dark Leaf Respiration (kgC/indiv/s) - - ! GPP IN THIS SUBROUTINE IS A RATE. THE CALLING ARGUMENT IS GPP_TSTEP. AFTER THIS - ! CALL THE RATE WILL BE MULTIPLIED BY THE INTERVAL TO GIVE THE INTEGRATED QUANT. - - - ! Locals - real(r8) :: tree_area - real(r8) :: laifrac - - ! Convert from umolC/m2leaf/s to umolC/indiv/s ( x canopy area x 1m2 leaf area). - tree_area = c_area/nplant - - ! The routine is only called if there are leaves. If there are leaves, - ! there is at least 1 layer - - laifrac = (treelai+treesai)-dble(nv-1)*dinc_ed - - ! Canopy Conductance - gscan = 1.0_r8/(rs_llz(nv)+rb)*laifrac*tree_area - - ! GPP - gpp = psn_llz(nv) * elai_llz(nv) * laifrac * tree_area - - ! Dark respiration - rdark = lmr_llz(nv) * elai_llz(nv) * laifrac * tree_area - - ! If there is more than one layer, add the sum over the others - if ( nv>1 ) then - gpp = gpp + sum(psn_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area - rdark = rdark + sum(lmr_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area - gscan = gscan + sum((1.0_r8/(rs_llz(1:nv-1) + rb ))) * tree_area - end if - - ! Convert dark respiration and GPP from umol/plant/s to kgC/plant/s - ! Note that + end subroutine LeafLayerPhotosynthesis - rdark = rdark * umolC_to_kgC - gpp = gpp * umolC_to_kgC - - - if ( DEBUG ) then - write(fates_log(),*) 'EDPhoto 816 ', gpp - write(fates_log(),*) 'EDPhoto 817 ', psn_llz(1:nv) - write(fates_log(),*) 'EDPhoto 820 ', nv - write(fates_log(),*) 'EDPhoto 821 ', elai_llz(1:nv) - write(fates_log(),*) 'EDPhoto 843 ', rdark - write(fates_log(),*) 'EDPhoto 871 ', laifrac - write(fates_log(),*) 'EDPhoto 872 ', tree_area - write(fates_log(),*) 'EDPhoto 873 ', nv - endif + ! ===================================================================================== + + subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv + psn_llz, & ! in %psn_z(cl,ft,1:currentCohort%nv) + lmr_llz, & ! in lmr_z(cl,ft,1:currentCohort%nv) + rs_llz, & ! in rs_z(cl,ft,1:currentCohort%nv) + anet_av_llz, & ! in anet_av_z(cl,ft,1:currentCohort%nv) + elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) + c_area, & ! in currentCohort%c_area + nplant, & ! in currentCohort%n + treelai, & ! in currentCohort%treelai + treesai, & ! in currentCohort%treesai + rb, & ! in bc_in(s)%rb_pa(ifp) + gscan, & ! out currentCohort%gscan + gpp, & ! out currentCohort%gpp_tstep + rdark) ! out currentCohort%rdark - return -end subroutine ScaleLeafLayerFluxToCohort + ! ------------------------------------------------------------------------------------ + ! This subroutine effectively integrates leaf carbon fluxes over the + ! leaf layers to give cohort totals. + ! Some arguments have the suffix "_llz". This indicates that the vector + ! is stratefied in the leaf-layer (ll) dimension, and is a portion of the calling + ! array which has the "_z" tag, thus "llz". + ! ------------------------------------------------------------------------------------ + + use FatesConstantsMod, only : umolC_to_kgC + use EDTypesMod, only : dinc_ed + + ! Arguments + integer, intent(in) :: nv ! number of active leaf layers + real(r8), intent(in) :: psn_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: lmr_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: rs_llz(nv) ! s/m + real(r8), intent(in) :: anet_av_llz(nv) ! umolC/m2leaf/s + real(r8), intent(in) :: elai_llz(nv) ! exposed LAI per layer + real(r8), intent(in) :: c_area ! crown area m2/m2 + real(r8), intent(in) :: nplant ! indiv/m2 + real(r8), intent(in) :: treelai ! m2/m2 + real(r8), intent(in) :: treesai ! m2/m2 + real(r8), intent(in) :: rb ! boundary layer resistance (s/m) + + real(r8), intent(out) :: gscan ! Canopy conductance of the cohort m/s + real(r8), intent(out) :: gpp ! GPP (kgC/indiv/s) + real(r8), intent(out) :: rdark ! Dark Leaf Respiration (kgC/indiv/s) + + ! GPP IN THIS SUBROUTINE IS A RATE. THE CALLING ARGUMENT IS GPP_TSTEP. AFTER THIS + ! CALL THE RATE WILL BE MULTIPLIED BY THE INTERVAL TO GIVE THE INTEGRATED QUANT. + + ! Locals + real(r8) :: tree_area + real(r8) :: laifrac + + ! Convert from umolC/m2leaf/s to umolC/indiv/s ( x canopy area x 1m2 leaf area). + tree_area = c_area/nplant + + ! The routine is only called if there are leaves. If there are leaves, + ! there is at least 1 layer + + laifrac = (treelai+treesai)-dble(nv-1)*dinc_ed + + ! Canopy Conductance + gscan = 1.0_r8/(rs_llz(nv)+rb)*laifrac*tree_area + + ! GPP + gpp = psn_llz(nv) * elai_llz(nv) * laifrac * tree_area + + ! Dark respiration + rdark = lmr_llz(nv) * elai_llz(nv) * laifrac * tree_area + + ! If there is more than one layer, add the sum over the others + if ( nv>1 ) then + gpp = gpp + sum(psn_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area + rdark = rdark + sum(lmr_llz(1:nv-1) * elai_llz(1:nv-1)) * tree_area + gscan = gscan + sum((1.0_r8/(rs_llz(1:nv-1) + rb ))) * tree_area + end if + + ! Convert dark respiration and GPP from umol/plant/s to kgC/plant/s + + rdark = rdark * umolC_to_kgC + gpp = gpp * umolC_to_kgC + + if ( DEBUG ) then + write(fates_log(),*) 'EDPhoto 816 ', gpp + write(fates_log(),*) 'EDPhoto 817 ', psn_llz(1:nv) + write(fates_log(),*) 'EDPhoto 820 ', nv + write(fates_log(),*) 'EDPhoto 821 ', elai_llz(1:nv) + write(fates_log(),*) 'EDPhoto 843 ', rdark + write(fates_log(),*) 'EDPhoto 871 ', laifrac + write(fates_log(),*) 'EDPhoto 872 ', tree_area + write(fates_log(),*) 'EDPhoto 873 ', nv + endif + + return + end subroutine ScaleLeafLayerFluxToCohort -! ======================================================================================= + ! ===================================================================================== -function ft1_f(tl, ha) result(ans) - ! - !!DESCRIPTION: + function ft1_f(tl, ha) result(ans) + ! + !!DESCRIPTION: ! photosynthesis temperature response ! ! !REVISION HISTORY @@ -1257,7 +1284,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) type(ed_cohort_type), pointer :: currentCohort ! Locals - integer :: CL ! Canopy Layer Index + integer :: cl ! Canopy Layer Index integer :: ft ! Function Type Index integer :: iv ! index of the exposed leaf layer for each canopy layer and pft @@ -1272,7 +1299,7 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) do while(associated(currentCohort)) currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) + max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) currentCohort => currentCohort%shorter @@ -1282,16 +1309,16 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) currentPatch%nrad = currentPatch%ncan ! Now loop through and identify which layer and pft combo has scattering elements - do CL = 1,cp_nclmax + do cl = 1,cp_nclmax do ft = 1,numpft_ed - currentPatch%present(CL,ft) = 0 - do iv = 1, currentPatch%nrad(CL,ft); - if(currentPatch%canopy_area_profile(CL,ft,iv) > 0._r8)then - currentPatch%present(CL,ft) = 1 + currentPatch%present(cl,ft) = 0 + do iv = 1, currentPatch%nrad(cl,ft); + if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then + currentPatch%present(cl,ft) = 1 end if end do !iv enddo !ft - enddo !CL + enddo !cl return end subroutine UpdateCanopyNCanNRadPresent @@ -1385,7 +1412,9 @@ subroutine GetCanopyGasParameters(can_press, & co2_cpoint = 1.0_r8 end if - ! THESE HARD CODED CONVERSIONS NEED TO BE CALLED FROM GLOBAL CONSTANTS (RGK 10-13-2016) + ! THESE HARD CODED CONVERSIONS NEED TO BE CALLED FROM GLOBAL CONSTANTS + ! (RGK 10-13-201). THE MEANING OF CF IS UNCLEAR, BUT THIS APPEARS TO BE A MOLAR CONVERSION + cf = can_press/(rgas*1.e-3_r8 * air_tempk )*1.e06_r8 gb_mol = (1._r8/ rb) * cf @@ -1410,22 +1439,22 @@ subroutine LeafLayerMaintenanceRespiration(lmr25top_ft, & use pftconMod , only : pftcon ! Arguments - real(r8), intent(in) :: lmr25top_ft ! canopy top leaf maint resp rate at 25C - ! for this pft (umol CO2/m**2/s) - integer, intent(in) :: ft ! (plant) Functional Type Index - real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile - real(r8), intent(in) :: veg_tempk ! vegetation temperature - real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) + real(r8), intent(in) :: lmr25top_ft ! canopy top leaf maint resp rate at 25C + ! for this pft (umol CO2/m**2/s) + integer, intent(in) :: ft ! (plant) Functional Type Index + real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile + real(r8), intent(in) :: veg_tempk ! vegetation temperature + real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) ! Locals - real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - + real(r8) :: lmr25 ! leaf layer: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) ! Parameter real(r8), parameter :: lmrha = 46390._r8 ! activation energy for lmr (J/mol) real(r8), parameter :: lmrhd = 150650._r8 ! deactivation energy for lmr (J/mol) real(r8), parameter :: lmrse = 490._r8 ! entropy term for lmr (J/mol/K) - real(r8), parameter :: lmrc = 1.15912391_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8), parameter :: lmrc = 1.15912391_r8 ! scaling factor for high + ! temperature inhibition (25 C = 1.0) ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s ! ---------------------------------------------------------------------------------- @@ -1442,7 +1471,7 @@ subroutine LeafLayerMaintenanceRespiration(lmr25top_ft, & ! Any hydrodynamic limitations could go here, currently none ! lmr = lmr * (nothing) - end subroutine LeafLayerMaintenanceRespiration + end subroutine LeafLayerMaintenanceRespiration ! ==================================================================================== @@ -1495,16 +1524,22 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & real(r8), intent(in) :: btran ! transpiration wetness factor (0 to 1) real(r8), intent(out) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8), intent(out) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8), intent(out) :: tpu ! triose phosphate utilization rate (umol CO2/m**2/s) + real(r8), intent(out) :: jmax ! maximum electron transport rate + ! (umol electrons/m**2/s) + real(r8), intent(out) :: tpu ! triose phosphate utilization rate + ! (umol CO2/m**2/s) real(r8), intent(out) :: co2_rcurve_islope ! initial slope of CO2 response curve (C4 plants) ! Locals ! ------------------------------------------------------------------------------- - real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C (umol CO2/m**2/s) - real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C (umol electrons/m**2/s) - real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C (umol CO2/m**2/s) - real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve (C4 plants) at 25C + real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C + ! (umol CO2/m**2/s) + real(r8) :: jmax25 ! leaf layer: maximum electron transport rate at 25C + ! (umol electrons/m**2/s) + real(r8) :: tpu25 ! leaf layer: triose phosphate utilization rate at 25C + ! (umol CO2/m**2/s) + real(r8) :: co2_rcurve_islope25 ! leaf layer: Initial slope of CO2 response curve + ! (C4 plants) at 25C ! Parameters @@ -1518,9 +1553,12 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & real(r8), parameter :: vcmaxse = 485._r8 ! entropy term for vcmax (J/mol/K) real(r8), parameter :: jmaxse = 495._r8 ! entropy term for jmax (J/mol/K) real(r8), parameter :: tpuse = 490._r8 ! entropy term for tpu (J/mol/K) - real(r8), parameter :: vcmaxc = 1.1534040_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8), parameter :: jmaxc = 1.1657242_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) - real(r8), parameter :: tpuc = 1.1591239_r8 ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8), parameter :: vcmaxc = 1.1534040_r8 ! scaling factor for high + ! temperature inhibition (25 C = 1.0) + real(r8), parameter :: jmaxc = 1.1657242_r8 ! scaling factor for high + ! temperature inhibition (25 C = 1.0) + real(r8), parameter :: tpuc = 1.1591239_r8 ! scaling factor for high + ! temperature inhibition (25 C = 1.0) if ( parsun_lsl <= 0._r8) then ! night time vcmax = 0._r8 @@ -1543,15 +1581,14 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & vcmax = vcmax / (1._r8 + exp( 0.2_r8*((tfrz+15._r8)-veg_tempk ) )) vcmax = vcmax / (1._r8 + exp( 0.3_r8*(veg_tempk-(tfrz+40._r8)) )) end if - co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) !q10 response of product limited psn. + !q10 response of product limited psn. + co2_rcurve_islope = co2_rcurve_islope25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) end if ! Adjust for water limitations vcmax = vcmax * btran - + return - end subroutine LeafLayerBiophysicalRates - - + end subroutine LeafLayerBiophysicalRates end module FATESPhotosynthesisMod From bf40dfd57dd79e8b9653bd4b7cfc883b0d1abd3e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 5 Dec 2016 15:23:26 -0800 Subject: [PATCH 07/15] Some minor format refactors, mostly line length truncation and removing unused variables. --- biogeophys/FatesPhotosynthesisMod.F90 | 170 +++++++++++++++----------- 1 file changed, 96 insertions(+), 74 deletions(-) diff --git a/biogeophys/FatesPhotosynthesisMod.F90 b/biogeophys/FatesPhotosynthesisMod.F90 index 7a313a6581..51ec13f3d5 100644 --- a/biogeophys/FatesPhotosynthesisMod.F90 +++ b/biogeophys/FatesPhotosynthesisMod.F90 @@ -103,6 +103,8 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! If we are not using hydraulics, we calculate a unique solution for each ! level-pft-layer combination. Thus the following three arrays are statically ! allocated for the maximum space of the two cases (numCohortsPerPatch) + ! The "_z" suffix indicates these variables are discretized at the "leaf_layer" + ! scale. ! ----------------------------------------------------------------------------------- ! leaf maintenance (dark) respiration (umol CO2/m**2/s) Double check this @@ -118,47 +120,52 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) ! used already logical :: rate_mask_z(cp_nclmax,mxpft,cp_nlevcan) - - real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation (umol co2/m**2/s) - real(r8) :: jmax_z ! leaf layer maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: tpu_z ! leaf layer triose phosphate utilization rate (umol CO2/m**2/s) - real(r8) :: kp_z ! leaf layer initial slope of CO2 response curve (C4 plants) - real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) - real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) - real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) - real(r8) :: co2_cpoint ! CO2 compensation point (Pa) - real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) - - - real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) - real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient - real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation at 25C (umol CO2/m**2/s) - real(r8) :: jmax25top(mxpft) ! canopy top: maximum electron transport rate at 25C (umol electrons/m**2/s) - real(r8) :: tpu25top(mxpft) ! canopy top: triose phosphate utilization rate at 25C (umol CO2/m**2/s) - real(r8) :: lmr25top(mxpft) ! canopy top: leaf maintenance respiration rate at 25C (umol CO2/m**2/s) - real(r8) :: kp25top(mxpft) ! canopy top: initial slope of CO2 response curve (C4 plants) at 25C - - ! Other - integer :: cl,s,iv,j,ps,ft,ifp ! indices - integer :: nv ! number of leaf layers - integer :: NCL_p ! number of canopy layers in patch - real(r8) :: cf ! s m**2/umol -> s/m - real(r8) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) - real(r8) :: ceair ! vapor pressure of air, constrained (Pa) - real(r8) :: nscaler ! leaf nitrogen scaling coefficient - real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass - real(r8) :: laican ! canopy sum of lai_z - real(r8) :: vai ! leaf and steam area in ths layer. - real(r8) :: laifrac - real(r8) :: tcsoi ! Temperature response function for root respiration. - real(r8) :: tcwood ! Temperature response function for wood - real(r8) :: tree_area - real(r8) :: gs_cohort - real(r8) :: rscanopy - real(r8) :: elai - 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) + real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation + ! (umol co2/m**2/s) + real(r8) :: jmax_z ! leaf layer maximum electron transport rate + ! (umol electrons/m**2/s) + real(r8) :: tpu_z ! leaf layer triose phosphate utilization rate + ! (umol CO2/m**2/s) + real(r8) :: kp_z ! leaf layer initial slope of CO2 response + ! curve (C4 plants) + real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) + real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) + real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) + real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) + real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient + real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation + ! at 25C (umol CO2/m**2/s) + real(r8) :: jmax25top(mxpft) ! canopy top: maximum electron transport + ! rate at 25C (umol electrons/m**2/s) + real(r8) :: tpu25top(mxpft) ! canopy top: triose phosphate utilization rate + ! at 25C (umol CO2/m**2/s) + real(r8) :: lmr25top(mxpft) ! canopy top: leaf maintenance respiration rate + ! at 25C (umol CO2/m**2/s) + real(r8) :: kp25top(mxpft) ! canopy top: initial slope of CO2 response curve + ! (C4 plants) at 25C + real(r8) :: cf ! s m**2/umol -> s/m + real(r8) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) + real(r8) :: ceair ! vapor pressure of air, constrained (Pa) + real(r8) :: nscaler ! leaf nitrogen scaling coefficient + real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass + real(r8) :: laican ! canopy sum of lai_z + real(r8) :: vai ! leaf and steam area in ths layer. + real(r8) :: tcsoi ! Temperature response function for root respiration. + real(r8) :: tcwood ! Temperature response function for wood + real(r8) :: rscanopy ! Canopy resistance [s/m] + real(r8) :: elai ! exposed LAI (patch scale) + 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) + + + integer :: cl,s,iv,j,ps,ft,ifp ! indices + integer :: nv ! number of leaf layers + integer :: NCL_p ! number of canopy layers in patch ! Parameters ! ----------------------------------------------------------------------- @@ -188,13 +195,15 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) associate( & - c3psn => pftcon%c3psn , & - slatop => pftcon%slatop , & ! specific leaf area at top of canopy, projected area basis [m^2/gC] - flnr => pftcon%flnr , & ! fraction of leaf N in the Rubisco enzyme (gN Rubisco / gN leaf) - 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) ! slope of BB relationship + c3psn => pftcon%c3psn , & + slatop => pftcon%slatop , & ! specific leaf area at top of canopy, + ! projected area basis [m^2/gC] + flnr => pftcon%flnr , & ! fraction of leaf N in the Rubisco + ! enzyme (gN Rubisco / gN leaf) + 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) ! slope of BB relationship q10 => EDParamsShareInst%Q10) do s = 1,nsites @@ -750,7 +759,6 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in real(r8), intent(out) :: anet_av_out ! net leaf photosynthesis (umol CO2/m**2/s) ! averaged over sun and shade leaves. - ! Locals ! ------------------------------------------------------------------------ integer :: pp_type ! Index for the different photosynthetic pathways C3,C4 @@ -772,7 +780,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in real(r8) :: gs_mol_err ! gs_mol for error check real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) - real(r8) :: ap ! product-limited (C3) or CO2-limited (C4) gross photosynthesis (umol CO2/m**2/s) + real(r8) :: ap ! product-limited (C3) or CO2-limited + ! (C4) gross photosynthesis (umol CO2/m**2/s) real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) ! Parameters @@ -826,14 +835,15 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. - rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. + rstoma_out = 0._r8 ! 1/rs is accumulated across sun and shaded leaves. anet_av_out = 0._r8 gstoma = 0._r8 do sunsha = 1,2 - ! Electron transport rate for C3 plants. Convert par from W/m2 to umol photons/m**2/s - ! using the factor 4.6 - ! Convert from units of par absorbed per unit ground area to par absorbed per unit leaf area. + ! Electron transport rate for C3 plants. + ! Convert par from W/m2 to umol photons/m**2/s using the factor 4.6 + ! Convert from units of par absorbed per unit ground area to par + ! absorbed per unit leaf area. if(sunsha == 1)then !sunlit if(( laisun_lsl * canopy_area_lsl) > 0.0000000001_r8)then @@ -885,7 +895,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in (co2_intra_c+mm_kco2 * (1._r8+can_o2_ppress / mm_ko2 )) ! C3: RuBP-limited photosynthesis - aj = je * max(co2_intra_c-co2_cpoint, 0._r8) / (4._r8*co2_intra_c+8._r8*co2_cpoint) + aj = je * max(co2_intra_c-co2_cpoint, 0._r8) / & + (4._r8*co2_intra_c+8._r8*co2_cpoint) ! C3: Product-limited photosynthesis ap = 3._r8 * tpu @@ -897,7 +908,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! C4: RuBP-limited photosynthesis if(sunsha == 1)then !sunlit - if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then !guard against /0's in the night. + !guard against /0's in the night. + if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then aj = quant_eff(pp_type) * parsun_lsl * 4.6_r8 !convert from per cohort to per m2 of leaf) aj = aj / (laisun_lsl * canopy_area_lsl) @@ -939,8 +951,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in leaf_co2_ppress = can_co2_ppress- 1.4_r8/gb_mol * anet * can_press leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) aquad = leaf_co2_ppress - bquad = leaf_co2_ppress*(gb_mol - bbb) - bb_slope(ft) * anet * can_press - cquad = -gb_mol*(leaf_co2_ppress*bbb + bb_slope(ft)*anet*can_press * ceair/ veg_esat ) + bquad = leaf_co2_ppress*(gb_mol - bbb) - bb_slope(ft) * anet * can_press + cquad = -gb_mol*(leaf_co2_ppress*bbb + & + bb_slope(ft)*anet*can_press * ceair/ veg_esat ) call quadratic_f (aquad, bquad, cquad, r1, r2) gs_mol = max(r1,r2) @@ -949,11 +962,13 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in co2_intra_c = can_co2_ppress - anet * can_press * & (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) - ! Check for co2_intra_c convergence. Delta co2_intra_c/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 + ! Check for co2_intra_c convergence. Delta co2_intra_c/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(co2_intra_c-co2_intra_c_old)/can_press*1.e06_r8 <= 2.e-06_r8) .or. niter == 5) then + if ((abs(co2_intra_c-co2_intra_c_old)/can_press*1.e06_r8 <= 2.e-06_r8) & + .or. niter == 5) then loop_continue = .false. end if end do !iteration loop @@ -963,10 +978,12 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in gs_mol = bbb end if - ! Final estimates for leaf_co2_ppress and co2_intra_c (needed for early exit of co2_intra_c iteration when an < 0) + ! Final estimates for leaf_co2_ppress and co2_intra_c + ! (needed for early exit of co2_intra_c iteration when an < 0) leaf_co2_ppress = can_co2_ppress - 1.4_r8/gb_mol * anet * can_press leaf_co2_ppress = max(leaf_co2_ppress,1.e-06_r8) - co2_intra_c = can_co2_ppress - anet * can_press * (1.4_r8*gs_mol+1.6_r8*gb_mol) / (gb_mol*gs_mol) + co2_intra_c = can_co2_ppress - anet * can_press * & + (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 @@ -975,7 +992,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross if ( DEBUG ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl - !accumulate total photosynthesis umol/m2 ground/s-1. weight per unit sun and sha leaves. + ! Accumulate total photosynthesis umol/m2 ground/s-1. + ! weight per unit sun and sha leaves. if(sunsha == 1)then !sunlit psn_out = psn_out + agross * f_sun_lsl anet_av_out = anet_av_out + anet * f_sun_lsl @@ -1012,7 +1030,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in !average leaf-level stomatal resistance rate over sun and shade leaves... rstoma_out = 1._r8/gstoma - else !No leaf area. This layer is present only because of stems. (leaves are off, or have reduced to 0 + else + !No leaf area. This layer is present only because of stems. + ! (leaves are off, or have reduced to 0) psn_out = 0._r8 rstoma_out = min(rsmax0, 1._r8/bbb * cf) @@ -1163,10 +1183,10 @@ function fth_f(tl,hd,se,scaleFactor) result(ans) ! ! !ARGUMENTS: - real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) - real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temperature function (J/mol) - real(r8), intent(in) :: se ! entropy term in photosynthesis temperature function (J/mol/K) - real(r8), intent(in) :: scaleFactor ! scaling factor for high temperature inhibition (25 C = 1.0) + real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temp function (K) + real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) + real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) + real(r8), intent(in) :: scaleFactor ! scaling factor for high temp inhibition (25 C = 1.0) ! ! !LOCAL VARIABLES: real(r8) :: ans @@ -1195,8 +1215,8 @@ function fth25_f(hd,se)result(ans) ! ! !ARGUMENTS: - real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temperature function (J/mol) - real(r8), intent(in) :: se ! entropy term in photosynthesis temperature function (J/mol/K) + real(r8), intent(in) :: hd ! deactivation energy in photosynthesis temp function (J/mol) + real(r8), intent(in) :: se ! entropy term in photosynthesis temp function (J/mol/K) ! ! !LOCAL VARIABLES: real(r8) :: ans @@ -1294,12 +1314,14 @@ subroutine UpdateCanopyNCanNRadPresent(currentPatch) ! --------------------------------------------------------------------------------- currentPatch%ncan(:,:) = 0 - !redo the canopy structure algorithm to get round a bug that is happening for site 125, FT13. + ! redo the canopy structure algorithm to get round a + ! bug that is happening for site 125, FT13. currentCohort => currentPatch%tallest do while(associated(currentCohort)) currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft),currentCohort%NV) + max(currentPatch%ncan(currentCohort%canopy_layer,currentCohort%pft), & + currentCohort%NV) currentCohort => currentCohort%shorter From 67cc470ea25a2ed8da617b6903d8b75f87c9ca98 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 6 Dec 2016 11:22:23 -0800 Subject: [PATCH 08/15] More minor syntactical changes to photosynthesis, some comments. --- ...od.F90 => FatesPlantRespPhotosynthMod.F90} | 21 ++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) rename biogeophys/{FatesPhotosynthesisMod.F90 => FatesPlantRespPhotosynthMod.F90} (98%) diff --git a/biogeophys/FatesPhotosynthesisMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 similarity index 98% rename from biogeophys/FatesPhotosynthesisMod.F90 rename to biogeophys/FatesPlantRespPhotosynthMod.F90 index 51ec13f3d5..a2b615a1f0 100644 --- a/biogeophys/FatesPhotosynthesisMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -1,8 +1,8 @@ -module FATESPhotosynthesisMod +module FATESPlantRespPhotosynthMod !------------------------------------------------------------------------------------- ! !DESCRIPTION: - ! Calculates the photosynthetic fluxes for the FATES model + ! Calculates the plant respiration and photosynthetic fluxes for the FATES model ! This code is similar to and was originally based off of the 'photosynthesis' ! subroutine in the CLM model. ! @@ -591,8 +591,12 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) 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 + + + ! ------------------------------------------------------------------ + ! Part IX: Perform some unit conversions (rate to integrated) and + ! calcualate some fluxes that are sums and nets of the base fluxes + ! ------------------------------------------------------------------ if ( DEBUG ) write(fates_log(),*) 'EDPhoto 904 ', currentCohort%resp_m if ( DEBUG ) write(fates_log(),*) 'EDPhoto 905 ', currentCohort%rdark @@ -630,6 +634,13 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) currentCohort%npp_tstep = currentCohort%gpp_tstep - & currentCohort%resp_tstep ! kgC/indiv/ts + + + ! psncanopy (gpp) and lmrcanopy (dark resp) are not used + ! by the host model right now. Once upon a time they were diagnostics. + ! Now we have our own diagnostics for GPP and LMR, so this step + ! is not really needed. + ! -------------------------------------------------------------------- bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + & currentCohort%gpp_tstep bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + & @@ -1613,4 +1624,4 @@ subroutine LeafLayerBiophysicalRates( parsun_lsl, & return end subroutine LeafLayerBiophysicalRates -end module FATESPhotosynthesisMod + end module FATESPlantRespPhotosynthMod From 98acd1ca04843cc4396e7eb39c04aa491cc7abbc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 6 Dec 2016 15:27:29 -0800 Subject: [PATCH 09/15] Changed the name of the photosynthesis module to reflect that it is fates, not ED, and that it handles respiration as well as photosynthesis. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index a2b615a1f0..362aea48e7 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -29,9 +29,7 @@ module FATESPlantRespPhotosynthMod implicit none private - - ! PUBLIC MEMBER FUNCTIONS: - public :: FATESPhotosynthesis !ED specific photosynthesis routine + public :: FatesPlantRespPhotosynthDrive ! Called by the HLM-Fates interface character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -46,7 +44,7 @@ module FATESPlantRespPhotosynthMod !-------------------------------------------------------------------------------------- - subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) + subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! ----------------------------------------------------------------------------------- ! !DESCRIPTION: @@ -681,7 +679,7 @@ subroutine FatesPhotosynthesis (nsites, sites,bc_in,bc_out,dtime) end do !site loop end associate - end subroutine FatesPhotosynthesis + end subroutine FatesPlantRespPhotosynthDrive ! ======================================================================================= From 386ded65d14f66a9c8130d83298782ed2f90934e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 9 Dec 2016 14:52:19 -0800 Subject: [PATCH 10/15] Fixed index order on the rate_mask_z array. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 362aea48e7..abe0c7deb6 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -375,7 +375,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! not been done yet. ! ------------------------------------------------------------ - if ( .not.rate_mask_z(ft,cl,iv) .or. use_fates_plant_hydro ) then + if ( .not.rate_mask_z(cl,ft,iv) .or. use_fates_plant_hydro ) then if (use_fates_plant_hydro) then write(fates_log(),*) 'use_fates_plant_hydro in EDTypes' @@ -466,7 +466,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) rs_z(cl,ft,iv), & ! out anet_av_z(cl,ft,iv)) ! out - rate_mask_z(ft,cl,iv) = .true. + rate_mask_z(cl,ft,iv) = .true. end if end do From 7a1c5f64ab6fc6bca947f756ccddb693dafa2d71 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 19 Dec 2016 14:54:06 -0800 Subject: [PATCH 11/15] Removed some unnecessary boundary conditions (LMR,PSN,GC). Reduced the scope of what is flushed in the photosynthesis mask to speed things up. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 134 +++++++-------------- main/EDTypesMod.F90 | 8 +- main/FatesInterfaceMod.F90 | 9 -- 3 files changed, 53 insertions(+), 98 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 362aea48e7..efb912b345 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -79,6 +79,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) use FatesConstantsMod, only : umol_per_mmol use FatesConstantsMod, only : rgas => rgas_J_K_kmol use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use FatesParameterDerivedMod, only : param_derived + ! ARGUMENTS: ! ----------------------------------------------------------------------------------- @@ -104,7 +106,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! The "_z" suffix indicates these variables are discretized at the "leaf_layer" ! scale. ! ----------------------------------------------------------------------------------- - + ! leaf maintenance (dark) respiration (umol CO2/m**2/s) Double check this real(r8) :: lmr_z(cp_nclmax,mxpft,cp_nlevcan) @@ -126,23 +128,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! (umol CO2/m**2/s) real(r8) :: kp_z ! leaf layer initial slope of CO2 response ! curve (C4 plants) - real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) + real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) real(r8) :: co2_cpoint ! CO2 compensation point (Pa) real(r8) :: btran_eff ! effective transpiration wetness factor (0 to 1) real(r8) :: bbb ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s) real(r8) :: kn(mxpft) ! leaf nitrogen decay coefficient - real(r8) :: vcmax25top(mxpft) ! canopy top: maximum rate of carboxylation - ! at 25C (umol CO2/m**2/s) - real(r8) :: jmax25top(mxpft) ! canopy top: maximum electron transport - ! rate at 25C (umol electrons/m**2/s) - real(r8) :: tpu25top(mxpft) ! canopy top: triose phosphate utilization rate - ! at 25C (umol CO2/m**2/s) - real(r8) :: lmr25top(mxpft) ! canopy top: leaf maintenance respiration rate - ! at 25C (umol CO2/m**2/s) - real(r8) :: kp25top(mxpft) ! canopy top: initial slope of CO2 response curve - ! (C4 plants) at 25C real(r8) :: cf ! s m**2/umol -> s/m real(r8) :: gb_mol ! leaf boundary layer conductance (umol H2O/m**2/s) real(r8) :: ceair ! vapor pressure of air, constrained (Pa) @@ -159,7 +151,14 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) 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) - + real(r8) :: gccanopy_pa ! Patch level canopy stomatal conductance [mmol m-2 s-1] + + ! ----------------------------------------------------------------------------------- + ! Keeping these two definitions in case they need to be added later + ! + ! ----------------------------------------------------------------------------------- + !real(r8) :: psncanopy_pa ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s) + !real(r8) :: lmrcanopy_pa ! patch sunlit leaf maintenance respiration rate (umol CO2/m**2/s) integer :: cl,s,iv,j,ps,ft,ifp ! indices integer :: nv ! number of leaf layers @@ -202,7 +201,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) fnitr => pftcon%fnitr , & ! foliage nitrogen limitation factor (-) leafcn => pftcon%leafcn , & ! leaf C:N (gC/gN) frootcn => pftcon%frootcn, & ! froot C:N (gc/gN) ! slope of BB relationship - q10 => EDParamsShareInst%Q10) + q10 => EDParamsShareInst%Q10 ) + do s = 1,nsites @@ -219,11 +219,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! Part I. Zero output boundary conditions ! --------------------------------------------------------------------------- - bc_out(s)%psncanopy_pa(ifp) = 0._r8 - bc_out(s)%lmrcanopy_pa(ifp) = 0._r8 bc_out(s)%rssun_pa(ifp) = 0._r8 bc_out(s)%rssha_pa(ifp) = 0._r8 - bc_out(s)%gccanopy_pa(ifp) = 0._r8 + + gccanopy_pa = 0._r8 + + !psncanopy_pa = 0._r8 + !lmrcanopy_pa = 0._r8 ! Part II. Filter out patches ! Patch level filter flag for photosynthesis calculations @@ -270,32 +272,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! ------------------------------------------------------------------------ do ft = 1,numpft_ed - - - ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) - lnc = 1._r8 / (slatop(ft) * leafcn(ft)) - - ! at the moment in ED we assume that there is no active N cycle. - ! This should change, of course. FIX(RF,032414) Sep2011. - ! fudge - shortcut using fnitr as a proxy for vcmax... - vcmax25top(ft) = fnitr(ft) - - ! Parameters derived from vcmax25top. - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of - ! Experimental Botany 44:907-920. 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)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) - - jmax25top(ft) = 1.67_r8 * vcmax25top(ft) - tpu25top(ft) = 0.167_r8 * vcmax25top(ft) - kp25top(ft) = 20000._r8 * vcmax25top(ft) - ! Nitrogen scaling factor. ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al ! (2010) Biogeosciences, 7, 1833-1859 @@ -304,20 +281,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if (bc_in(s)%dayl_factor_pa(ifp) == 0._r8) then kn(ft) = 0._r8 else - kn(ft) = exp(0.00963_r8 * vcmax25top(ft) - 2.43_r8) + kn(ft) = exp(0.00963_r8 * param_derived%vcmax25top(ft) - 2.43_r8) end if - - ! Leaf maintenance respiration to match the base rate used in CN - ! but with the new temperature functions for C3 and C4 plants. - ! - ! - ! 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 - ! - ! 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 / (umolC_to_kgC * g_per_kg) end do !ft @@ -339,8 +304,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! With plant hydraulics, we must realize that photosynthesis and ! respiration will be different for leaves of each cohort in the leaf ! layers, as they will have there own hydraulic limitations. + ! NOTE: Only need to flush mask on the number of used pfts, not the whole + ! scratch space. ! ------------------------------------------------------------------------ - rate_mask_z(:,:,:) = .false. + rate_mask_z(:,1:numpft_ed,:) = .false. if(currentPatch%countcohorts > 0.0)then ! Ignore empty patches @@ -375,7 +342,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! not been done yet. ! ------------------------------------------------------------ - if ( .not.rate_mask_z(ft,cl,iv) .or. use_fates_plant_hydro ) then + if ( .not.rate_mask_z(cl,ft,iv) .or. use_fates_plant_hydro ) then if (use_fates_plant_hydro) then write(fates_log(),*) 'use_fates_plant_hydro in EDTypes' @@ -404,7 +371,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) nscaler = exp(-kn(ft) * laican) ! Part VII: Calculate dark respiration (leaf maintenance) for this layer - call LeafLayerMaintenanceRespiration( lmr25top(ft), & ! in + call LeafLayerMaintenanceRespiration( param_derived%lmr25top(ft),& ! in nscaler, & ! in ft, & ! in bc_in(s)%t_veg_pa(ifp), & ! in @@ -422,10 +389,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) call LeafLayerBiophysicalRates(currentPatch%ed_parsun_z(cl,ft,iv), & ! in ft, & ! in - vcmax25top(ft), & ! in - jmax25top(ft), & ! in - tpu25top(ft), & ! in - kp25top(ft), & ! in + param_derived%vcmax25top(ft), & ! in + param_derived%jmax25top(ft), & ! in + param_derived%tpu25top(ft), & ! in + param_derived%kp25top(ft), & ! in nscaler, & ! in bc_in(s)%t_veg_pa(ifp), & ! in btran_eff, & ! in @@ -466,7 +433,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) rs_z(cl,ft,iv), & ! out anet_av_z(cl,ft,iv)) ! out - rate_mask_z(ft,cl,iv) = .true. + rate_mask_z(cl,ft,iv) = .true. end if end do @@ -631,25 +598,15 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) currentCohort%resp_g ! kgC/indiv/ts currentCohort%npp_tstep = currentCohort%gpp_tstep - & currentCohort%resp_tstep ! kgC/indiv/ts + + ! Accumulate stomatal conductance over the patch + gccanopy_pa = gccanopy_pa + & + currentCohort%gscan * & + currentCohort%n /currentPatch%total_canopy_area - - - ! psncanopy (gpp) and lmrcanopy (dark resp) are not used - ! by the host model right now. Once upon a time they were diagnostics. - ! Now we have our own diagnostics for GPP and LMR, so this step - ! is not really needed. - ! -------------------------------------------------------------------- - bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) + & - currentCohort%gpp_tstep - bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) + & - currentCohort%resp_m - - ! accumulate cohort level canopy conductances over - ! whole area before dividing by total area - bc_out(s)%gccanopy_pa(ifp) = bc_out(s)%gccanopy_pa(ifp) + & - currentCohort%gscan * & - currentCohort%n /currentPatch%total_canopy_area - + !psncanopy_pa = psncanopy_pa + currentCohort%gpp_tstep + !lmrcanopy_pa = lmrcanopy_pa + currentCohort%resp_m + currentCohort => currentCohort%shorter enddo ! end cohort loop. @@ -658,18 +615,19 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) elai = calc_areaindex(currentPatch,'elai') - bc_out(s)%psncanopy_pa(ifp) = bc_out(s)%psncanopy_pa(ifp) / currentPatch%area - bc_out(s)%lmrcanopy_pa(ifp) = bc_out(s)%lmrcanopy_pa(ifp) / currentPatch%area - - if(bc_out(s)%gccanopy_pa(ifp) > 1._r8/rsmax0 .and. elai > 0.0_r8)then - rscanopy = (1.0_r8/bc_out(s)%gccanopy_pa(ifp))-bc_in(s)%rb_pa(ifp)/elai + !psncanopy_pa(ifp) = psncanopy_pa(ifp) / currentPatch%area + !lmrcanopy_pa(ifp) = lmrcanopy_pa(ifp) / currentPatch%area + + if(gccanopy_pa > 1._r8/rsmax0 .and. elai > 0.0_r8)then + rscanopy = (1.0_r8/gccanopy_pa)-bc_in(s)%rb_pa(ifp)/elai else rscanopy = rsmax0 end if + bc_out(s)%rssun_pa(ifp) = rscanopy bc_out(s)%rssha_pa(ifp) = rscanopy - !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 + + end if currentPatch => currentPatch%younger diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index d3e697e309..5896931f5c 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -614,7 +614,7 @@ function map_clmpatch_to_edpatch(site, clmpatch_number) result(edpatch_pointer) type(ed_patch_type), pointer :: edpatch_pointer !---------------------------------------------------------------------- - ! There is a one-to-one mapping between edpatches and clmpatches. To obtain + ! There is a one-to-one mapping between edpatches and clmpatches. To obtainFatesDerivedFromParameters ! this mapping - the following is computed elsewhere in the code base ! (1) what is the weight respective to the column of clmpatch? ! dynEDMod determines this via the following logic @@ -670,4 +670,10 @@ subroutine set_root_fraction( this , depth_gl) end subroutine set_root_fraction + + ! ===================================================================================== + + + + end module EDTypesMod diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 356951bcf0..752e35762f 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -191,9 +191,6 @@ module FatesInterfaceMod ! Shaded canopy resistance [s/m] real(r8), allocatable :: rssha_pa(:) - ! Canopy conductance [mmol m-2 s-1] - real(r8), allocatable :: gccanopy_pa(:) - ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s) real(r8), allocatable :: psncanopy_pa(:) @@ -406,9 +403,6 @@ subroutine allocate_bcout(bc_out) ! Photosynthesis allocate(bc_out%rssun_pa(numPatchesPerCol)) allocate(bc_out%rssha_pa(numPatchesPerCol)) - allocate(bc_out%gccanopy_pa(numPatchesPerCol)) - allocate(bc_out%lmrcanopy_pa(numPatchesPerCol)) - allocate(bc_out%psncanopy_pa(numPatchesPerCol)) ! Canopy Radiation allocate(bc_out%albd_parb(numPatchesPerCol,cp_numSWb)) @@ -481,9 +475,6 @@ subroutine zero_bcs(this,s) this%bc_out(s)%rssun_pa(:) = 0.0_r8 this%bc_out(s)%rssha_pa(:) = 0.0_r8 - this%bc_out(s)%gccanopy_pa(:) = 0.0_r8 - this%bc_out(s)%psncanopy_pa(:) = 0.0_r8 - this%bc_out(s)%lmrcanopy_pa(:) = 0.0_r8 this%bc_out(s)%albd_parb(:,:) = 0.0_r8 this%bc_out(s)%albi_parb(:,:) = 0.0_r8 From 533901c280f29158a07f9f4524d002a9ec211d61 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 19 Dec 2016 14:55:50 -0800 Subject: [PATCH 12/15] Correction of minor typo. --- main/EDTypesMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 5896931f5c..c70f5c079d 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -614,7 +614,7 @@ function map_clmpatch_to_edpatch(site, clmpatch_number) result(edpatch_pointer) type(ed_patch_type), pointer :: edpatch_pointer !---------------------------------------------------------------------- - ! There is a one-to-one mapping between edpatches and clmpatches. To obtainFatesDerivedFromParameters + ! There is a one-to-one mapping between edpatches and clmpatches. To obtain ! this mapping - the following is computed elsewhere in the code base ! (1) what is the weight respective to the column of clmpatch? ! dynEDMod determines this via the following logic From c7ed1ca332780856b6d01f865b7f3210101f0275 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 20 Dec 2016 14:56:48 -0800 Subject: [PATCH 13/15] Adding the parameter derived module to version control. --- main/FatesParameterDerivedMod.F90 | 116 ++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 main/FatesParameterDerivedMod.F90 diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 new file mode 100644 index 0000000000..41641d754e --- /dev/null +++ b/main/FatesParameterDerivedMod.F90 @@ -0,0 +1,116 @@ +module FatesParameterDerivedMod + + ! ------------------------------------------------------------------------------------- + ! This module contains all procedures types and settings for any quantities that are + ! statically derived from static model parameters. These are unchanging quantities + ! and are based off of simple relationships from parameters that the user can + ! vary. This should be called once, and early in the model initialization call + ! sequence immediately after FATES parameters are read in. + ! ------------------------------------------------------------------------------------- + + use FatesConstantsMod, only : r8 => fates_r8 + use FatesConstantsMod, only : umolC_to_kgC + use FatesConstantsMod, only : g_per_kg + + type param_derived_type + + real(r8), allocatable :: vcmax25top(:) ! canopy top: maximum rate of carboxylation + ! at 25C (umol CO2/m**2/s) + real(r8), allocatable :: jmax25top(:) ! canopy top: maximum electron transport + ! rate at 25C (umol electrons/m**2/s) + real(r8), allocatable :: tpu25top(:) ! canopy top: triose phosphate utilization + ! rate at 25C (umol CO2/m**2/s) + real(r8), allocatable :: kp25top(:) ! canopy top: initial slope of CO2 response + ! curve (C4 plants) at 25C + real(r8), allocatable :: lmr25top(:) ! canopy top: leaf maintenance respiration + ! rate at 25C (umol CO2/m**2/s) + contains + + procedure :: Init + procedure :: InitAllocate + + end type param_derived_type + + type(param_derived_type) :: param_derived + +contains + + subroutine InitAllocate(this,maxpft) + + class(param_derived_type), intent(inout) :: this + integer, intent(in) :: maxpft + + allocate(this%vcmax25top(maxpft)) + allocate(this%jmax25top(maxpft)) + allocate(this%tpu25top(maxpft)) + allocate(this%kp25top(maxpft)) + allocate(this%lmr25top(maxpft)) + + return + end subroutine InitAllocate + + ! ===================================================================================== + + subroutine Init(this,maxpft) + + use pftconMod , only: pftcon + + class(param_derived_type), intent(inout) :: this + integer, intent(in) :: maxpft + + ! local variables + integer :: ft ! pft index + real(r8) :: lnc ! leaf N concentration (gN leaf/m^2) + + associate( & + slatop => pftcon%slatop , & ! specific leaf area at top of canopy, + ! projected area basis [m^2/gC] + fnitr => pftcon%fnitr , & ! foliage nitrogen limitation factor (-) + leafcn => pftcon%leafcn ) ! leaf C:N (gC/gN) + + call this%InitAllocate(maxpft) + + do ft = 1,maxpft + + ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) + lnc = 1._r8 / (slatop(ft) * leafcn(ft)) + + ! at the moment in ED we assume that there is no active N cycle. + ! This should change, of course. FIX(RF,032414) Sep2011. + ! fudge - shortcut using fnitr as a proxy for vcmax... + this%vcmax25top(ft) = fnitr(ft) + + ! Parameters derived from vcmax25top. + ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 + ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of + ! Experimental Botany 44:907-920. 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)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) + + this%jmax25top(ft) = 1.67_r8 * this%vcmax25top(ft) + this%tpu25top(ft) = 0.167_r8 * this%vcmax25top(ft) + this%kp25top(ft) = 20000._r8 * this%vcmax25top(ft) + + ! Leaf maintenance respiration to match the base rate used in CN + ! but with the new temperature functions for C3 and C4 plants. + ! + ! + ! 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 + ! + ! Then scale this value at the top of the canopy for canopy depth + + this%lmr25top(ft) = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) + this%lmr25top(ft) = this%lmr25top(ft) * lnc / (umolC_to_kgC * g_per_kg) + + end do !ft + end associate + return + end subroutine Init + + +end module FatesParameterDerivedMod From 5908004e5d3283f9ab126c7be6cbfbe4eb3880bc Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 21 Dec 2016 11:40:23 -0800 Subject: [PATCH 14/15] Some minor optimizations in leaf layer photosynthesis; including the commenting out of debug calls, setting the initialization of a loop earlier in the code to remove a logic call, and adding an optional fast quadratic solver that simply removes error checking. The fast solve was left out, as this PR changegroup is not targetting optimization. --- biogeophys/FatesPlantRespPhotosynthMod.F90 | 88 ++++++++++++++++------ 1 file changed, 66 insertions(+), 22 deletions(-) diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index efb912b345..2b15e137be 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -346,9 +346,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if (use_fates_plant_hydro) then write(fates_log(),*) 'use_fates_plant_hydro in EDTypes' - write(fates_log(),*) 'has been set to true. You have inadvartently' + write(fates_log(),*) 'has been set to true. You have inadvertently' write(fates_log(),*) 'turned on a future feature that is not in the' - write(fates_log(),*) 'FATES model codeset yet. Please set this to' + write(fates_log(),*) 'FATES codeset yet. Please set this to' write(fates_log(),*) 'false and re-compile.' call endrun(msg=errMsg(sourcefile, __LINE__)) !! !! bbb = max (bbbopt(ps)*currentCohort%btran(iv), 1._r8) @@ -750,7 +750,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in real(r8) :: ap ! product-limited (C3) or CO2-limited ! (C4) gross photosynthesis (umol CO2/m**2/s) real(r8) :: ai ! intermediate co-limited photosynthesis (umol CO2/m**2/s) - real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + real(r8) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + real(r8) :: init_co2_intra_c ! First guess intracellular co2 specific to C path ! Parameters ! ------------------------------------------------------------------------ ! Fraction of light absorbed by non-photosynthetic pigments @@ -777,8 +778,10 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in if (nint(pftcon%c3psn(ft)) == 1) then! photosynthetic pathway: 0. = c4, 1. = c3 pp_type = 1 + init_co2_intra_c = init_a2l_co2_c3 * can_co2_ppress else pp_type = 2 + init_co2_intra_c = init_a2l_co2_c4 * can_co2_ppress end if ! Part III: Photosynthesis and Conductance @@ -792,13 +795,13 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in else ! day time (a little bit more complicated ...) - if ( DEBUG ) write(fates_log(),*) 'EDphot 594 ',laisun_lsl - if ( DEBUG ) write(fates_log(),*) 'EDphot 595 ',laisha_lsl +! if ( DEBUG ) write(fates_log(),*) 'EDphot 594 ',laisun_lsl +! if ( DEBUG ) write(fates_log(),*) 'EDphot 595 ',laisha_lsl !is there leaf area? - (NV can be larger than 0 with only stem area if deciduous) if ( laisun_lsl + laisha_lsl > 0._r8 ) then - if ( DEBUG ) write(fates_log(),*) '600 in laisun, laisha loop ' +! if ( DEBUG ) write(fates_log(),*) '600 in laisun, laisha loop ' !Loop aroun shaded and unshaded leaves psn_out = 0._r8 ! psn is accumulated across sun and shaded leaves. @@ -836,14 +839,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in call quadratic_f (aquad, bquad, cquad, r1, r2) je = min(r1,r2) - ! Iterative loop for ci beginning with initial guess - ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK) - - if (pp_type == 1)then - co2_intra_c = init_a2l_co2_c3 * can_co2_ppress - else - co2_intra_c = init_a2l_co2_c4 * can_co2_ppress - end if + ! Initialize intracellular co2 + co2_intra_c = init_co2_intra_c niter = 0 loop_continue = .true. @@ -855,7 +852,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in co2_intra_c_old = co2_intra_c ! Photosynthesis limitation rate calculations - if (pp_type == 1)then + if (pp_type == 1)then ! C3: Rubisco-limited photosynthesis ac = vcmax * max(co2_intra_c-co2_cpoint, 0._r8) / & @@ -875,7 +872,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! C4: RuBP-limited photosynthesis if(sunsha == 1)then !sunlit - !guard against /0's in the night. + !guard against /0's in the night. if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then aj = quant_eff(pp_type) * parsun_lsl * 4.6_r8 !convert from per cohort to per m2 of leaf) @@ -955,9 +952,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in ! 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(fates_log(),*) 'EDPhoto 737 ', psn_out - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 737 ', psn_out +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 738 ', agross +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 739 ', f_sun_lsl ! Accumulate total photosynthesis umol/m2 ground/s-1. ! weight per unit sun and sha leaves. @@ -972,9 +969,9 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in 1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl) end if - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', psn_out - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', agross - if ( DEBUG ) write(fates_log(),*) 'EDPhoto 760 ', f_sun_lsl +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 758 ', psn_out +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 759 ', agross +! if ( DEBUG ) write(fates_log(),*) 'EDPhoto 760 ', f_sun_lsl ! Make sure iterative solution is correct if (gs_mol < 0._r8) then @@ -1238,6 +1235,53 @@ subroutine quadratic_f (a, b, c, r1, r2) end if end subroutine quadratic_f + + ! ==================================================================================== + + subroutine quadratic_fast (a, b, c, r1, r2) + ! + ! !DESCRIPTION: + !==============================================================================! + !----------------- Solve quadratic equation for its two roots -----------------! + ! THIS METHOD SIMPLY REMOVES THE DIV0 CHECK AND ERROR REPORTING ! + !==============================================================================! + ! Solution from Press et al (1986) Numerical Recipes: The Art of Scientific + ! Computing (Cambridge University Press, Cambridge), pp. 145. + ! + ! !REVISION HISTORY: + ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson + ! 7/23/16: Copied over from CLM by Ryan Knox + ! + ! !USES: + ! + ! !ARGUMENTS: + real(r8), intent(in) :: a,b,c ! Terms for quadratic equation + real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation + ! + ! !LOCAL VARIABLES: + real(r8) :: q ! Temporary term for quadratic solution + !------------------------------------------------------------------------------ + + ! if (a == 0._r8) then + ! write (fates_log(),*) 'Quadratic solution error: a = ',a + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + ! end if + + if (b >= 0._r8) then + q = -0.5_r8 * (b + sqrt(b*b - 4._r8*a*c)) + else + q = -0.5_r8 * (b - sqrt(b*b - 4._r8*a*c)) + end if + + r1 = q / a + ! if (q /= 0._r8) then + r2 = c / q + ! else + ! r2 = 1.e36_r8 + ! end if + + end subroutine quadratic_fast + ! ==================================================================================== From 14bb4c870fc118a49bd93be5d8951702cb786a0a Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 9 Jan 2017 14:06:31 -0800 Subject: [PATCH 15/15] Small bug fix on previous merge conflict. Changed index ordering on temp photosynthesis arrays for performance. Added in a fix currently in another PR to help with testing (npp_acc in copy cohort). --- biogeochem/EDCohortDynamicsMod.F90 | 4 ++- biogeophys/FatesPlantRespPhotosynthMod.F90 | 41 ++++++++++++---------- main/FatesInterfaceMod.F90 | 4 +-- 3 files changed, 28 insertions(+), 21 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 3344cbb969..8d4f4bd553 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1037,13 +1037,15 @@ subroutine copy_cohort( currentCohort,copyc ) n%gpp_acc_hold = o%gpp_acc_hold n%gpp_acc = o%gpp_acc n%gpp_tstep = o%gpp_tstep + n%npp_acc_hold = o%npp_acc_hold n%npp_tstep = o%npp_tstep if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ia ',o%npp_acc if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ib ',o%resp_acc - n%npp_acc_hold = o%npp_acc_hold + n%npp_acc = o%npp_acc + n%resp_tstep = o%resp_tstep n%resp_acc = o%resp_acc n%resp_acc_hold = o%resp_acc_hold diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 2b15e137be..73f995df4d 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -105,20 +105,25 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! allocated for the maximum space of the two cases (numCohortsPerPatch) ! The "_z" suffix indicates these variables are discretized at the "leaf_layer" ! scale. + ! Note: For these temporary arrays, we have the leaf layer dimension first + ! and the canopy layer last. This order is chosen for efficiency. The arrays + ! such as leaf area that are bound to the patch structure DO NOT follow this order + ! as they are used in many other parts of the code with different looping, we + ! are not modifying its order now. ! ----------------------------------------------------------------------------------- ! leaf maintenance (dark) respiration (umol CO2/m**2/s) Double check this - real(r8) :: lmr_z(cp_nclmax,mxpft,cp_nlevcan) + real(r8) :: lmr_z(cp_nlevcan,mxpft,cp_nclmax) ! stomatal resistance s/m - real(r8) :: rs_z(cp_nclmax,mxpft,cp_nlevcan) + real(r8) :: rs_z(cp_nlevcan,mxpft,cp_nclmax) ! net leaf photosynthesis averaged over sun and shade leaves. (umol CO2/m**2/s) - real(r8) :: anet_av_z(cp_nclmax,mxpft,cp_nlevcan) + real(r8) :: anet_av_z(cp_nlevcan,mxpft,cp_nclmax) ! Mask used to determine which leaf-layer biophysical rates have been ! used already - logical :: rate_mask_z(cp_nclmax,mxpft,cp_nlevcan) + logical :: rate_mask_z(cp_nlevcan,mxpft,cp_nclmax) real(r8) :: vcmax_z ! leaf layer maximum rate of carboxylation ! (umol co2/m**2/s) @@ -342,7 +347,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! not been done yet. ! ------------------------------------------------------------ - if ( .not.rate_mask_z(cl,ft,iv) .or. use_fates_plant_hydro ) then + if ( .not.rate_mask_z(iv,ft,cl) .or. use_fates_plant_hydro ) then if (use_fates_plant_hydro) then write(fates_log(),*) 'use_fates_plant_hydro in EDTypes' @@ -375,7 +380,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) nscaler, & ! in ft, & ! in bc_in(s)%t_veg_pa(ifp), & ! in - lmr_z(cl,ft,iv)) ! out + lmr_z(iv,ft,cl)) ! out ! Part VII: Calculate (1) maximum rate of carboxylation (vcmax), ! (2) maximum electron transport rate, (3) triose phosphate @@ -428,12 +433,12 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) mm_kco2, & ! in mm_ko2, & ! in co2_cpoint, & ! in - lmr_z(cl,ft,iv), & ! in + lmr_z(iv,ft,cl), & ! in currentPatch%psn_z(cl,ft,iv), & ! out - rs_z(cl,ft,iv), & ! out - anet_av_z(cl,ft,iv)) ! out + rs_z(iv,ft,cl), & ! out + anet_av_z(iv,ft,cl)) ! out - rate_mask_z(cl,ft,iv) = .true. + rate_mask_z(iv,ft,cl) = .true. end if end do @@ -455,9 +460,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) nv = currentCohort%nv call ScaleLeafLayerFluxToCohort(nv, & !in currentPatch%psn_z(cl,ft,1:nv), & !in - lmr_z(cl,ft,1:nv), & !in - rs_z(cl,ft,1:nv), & !in - anet_av_z(cl,ft,1:nv), & !in + lmr_z(1:nv,ft,cl), & !in + rs_z(1:nv,ft,cl), & !in + anet_av_z(1:nv,ft,cl), & !in currentPatch%elai_profile(cl,ft,1:nv), & !in currentCohort%c_area, & !in currentCohort%n, & !in @@ -469,7 +474,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) currentCohort%rdark) !out ! Net Uptake does not need to be scaled, just transfer directly - currentCohort%ts_net_uptake(1:nv) = anet_av_z(cl,ft,1:nv) * umolC_to_kgC + currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC else @@ -1011,10 +1016,10 @@ end subroutine LeafLayerPhotosynthesis ! ===================================================================================== subroutine ScaleLeafLayerFluxToCohort(nv, & ! in currentCohort%nv - psn_llz, & ! in %psn_z(cl,ft,1:currentCohort%nv) - lmr_llz, & ! in lmr_z(cl,ft,1:currentCohort%nv) - rs_llz, & ! in rs_z(cl,ft,1:currentCohort%nv) - anet_av_llz, & ! in anet_av_z(cl,ft,1:currentCohort%nv) + psn_llz, & ! in %psn_z(1:currentCohort%nv,ft,cl) + lmr_llz, & ! in lmr_z(1:currentCohort%nv,ft,cl) + rs_llz, & ! in rs_z(1:currentCohort%nv,ft,cl) + anet_av_llz, & ! in anet_av_z(1:currentCohort%nv,ft,cl) elai_llz, & ! in %elai_profile(cl,ft,1:currentCohort%nv) c_area, & ! in currentCohort%c_area nplant, & ! in currentCohort%n diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 06e26b3771..139dfb7fac 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -404,8 +404,8 @@ subroutine allocate_bcout(bc_out) ! Photosynthesis - allocate(bc_out%rssun_pa(numPatchesPerCol)) - allocate(bc_out%rssha_pa(numPatchesPerCol)) + allocate(bc_out%rssun_pa(maxPatchesPerCol)) + allocate(bc_out%rssha_pa(maxPatchesPerCol)) ! Canopy Radiation allocate(bc_out%albd_parb(maxPatchesPerCol,cp_numSWb))