diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index d6ae9f52c8..98c39ffc22 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -11,6 +11,7 @@ module EDCanopyStructureMod use FatesConstantsMod , only : nearzero, area_error_1 use FatesConstantsMod , only : rsnbl_math_prec use FatesConstantsMod , only : nocomp_bareground + use FatesConstantsMod , only : fates_unset_r8 use FatesConstantsMod, only : i_term_mort_type_canlev use FatesGlobals , only : fates_log use EDPftvarcon , only : EDPftvarcon_inst @@ -23,12 +24,16 @@ module EDCanopyStructureMod use EDTypesMod , only : set_patchno use FatesAllometryMod , only : VegAreaLayer use FatesAllometryMod , only : CrownDepth + use FatesAllometryMod , only : bleaf + use FatesAllometryMod, only : storage_fraction_of_target + use LeafBiophysicsMod, only : LowstorageMainRespReduction use FatesPatchMod, only : fates_patch_type use FatesCohortMod, only : fates_cohort_type use EDParamsMod , only : nclmax use EDParamsMod , only : nlevleaf use EDParamsMod , only : GetNVegLayers use EDParamsMod , only : comp_excln_exp + use EDParamsMod , only : maxpft use EDtypesMod , only : AREA use EDLoggingMortalityMod , only : UpdateHarvestC use FatesGlobals , only : endrun => fates_endrun @@ -37,7 +42,8 @@ module EDCanopyStructureMod use FatesInterfaceTypesMod , only : hlm_use_cohort_age_tracking use FatesInterfaceTypesMod , only : hlm_use_sp use FatesInterfaceTypesMod , only : numpft - use FatesInterfaceTypesMod, only : bc_in_type + use FatesInterfaceTypesMod, only : bc_in_type + use FatesInterfaceTypesMod , only : fates_maxPatchesPerSite use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : leaf_organ @@ -51,6 +57,11 @@ module EDCanopyStructureMod use FatesTwoStreamUtilsMod, only : FatesConstructRadElements use FatesRadiationMemMod , only : twostr_solver use FatesRadiationMemMod , only : num_rad_stream_types + use LeafBiophysicsMod , only : UpdateSlowBiophysicalRates + use LeafBiophysicsMod, only : DecayCoeffVcmax + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use FatesInterfaceTypesMod, only : hlm_parteh_mode ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg @@ -349,7 +360,7 @@ subroutine canopy_structure( currentSite , bc_in ) write(fates_log(),*) 'and you think this number of canopy layers is reasonable.' call endrun(msg=errMsg(sourcefile, __LINE__)) else - currentPatch%NCL_p = z + currentPatch%ncl = z end if ! ------------------------------------------------------------------------------------------- @@ -689,11 +700,11 @@ subroutine PromoteOrDemote(site,patch,target_layer,phase,target_area) ! keep track of number and biomass promoted/demoted if( layer_co(ic)%pd_area > 0._r8 ) then - leaf_c = cohort%prt%GetState(leaf_organ,carbon12_element) - store_c = cohort%prt%GetState(store_organ,carbon12_element) - fnrt_c = cohort%prt%GetState(fnrt_organ,carbon12_element) - sapw_c = cohort%prt%GetState(sapw_organ,carbon12_element) - struct_c = cohort%prt%GetState(struct_organ,carbon12_element) + leaf_c = sum(cohort%prt%leaf_c(:)) + store_c = cohort%prt%store_c + fnrt_c = cohort%prt%fnrt_c + sapw_c = cohort%prt%sapw_c + struct_c = cohort%prt%struct_c if(phase==demotion_phase) then site%demotion_rate(cohort%size_class) = & @@ -787,7 +798,8 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDtypesMod , only : area use FatesConstantsMod , only : itrue - + use EDParamsMod , only : maxpatch_total + ! !ARGUMENTS integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) @@ -798,6 +810,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) type (fates_cohort_type) , pointer :: currentCohort integer :: s integer :: ft ! plant functional type + integer :: ift ! index of unique plant functional type array integer :: ifp ! the number of the vegetated patch (1,2,3). In SP mode bareground patch is 0 integer :: patchn ! identification number for each patch. real(r8) :: leaf_c ! leaf carbon [kg] @@ -805,7 +818,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] - + !---------------------------------------------------------------------- if ( debug ) then @@ -813,32 +826,44 @@ subroutine canopy_summarization( nsites, sites, bc_in ) endif do s = 1,nsites - + ! -------------------------------------------------------------------------------- ! Set the patch indices (this is usefull mostly for communicating with a host or ! driving model. Loops through all patches and sets cpatch%patchno to the integer ! order of oldest to youngest where the oldest is 1. ! -------------------------------------------------------------------------------- call set_patchno( sites(s) , .false., 0) - + currentPatch => sites(s)%oldest_patch - do while(associated(currentPatch)) - + !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 currentPatch%total_tree_area = 0.0_r8 + ! Tracking unique pft counts on this patch + ! used to help speed up photosynthesis calculations + currentPatch%unique_pfts(:) = 0 + currentPatch%upft_index(:) = 0 + currentPatch%nupft = 0 + !update cohort quantitie s currentCohort => currentPatch%shortest - do while(associated(currentCohort)) + co_loop1: do while(associated(currentCohort)) ft = currentCohort%pft - leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) - sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) - struct_c = currentCohort%prt%GetState(struct_organ, carbon12_element) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, carbon12_element) - store_c = currentCohort%prt%GetState(store_organ, carbon12_element) + + if( .not. any(currentPatch%unique_pfts(1:numpft) == ft) )then + currentPatch%nupft = currentPatch%nupft+1 + currentPatch%unique_pfts(currentPatch%nupft) = ft + currentPatch%upft_index(ft) = currentPatch%nupft + end if + + leaf_c = sum(currentCohort%prt%leaf_c(:)) + sapw_c = currentCohort%prt%sapw_c + struct_c = currentCohort%prt%struct_c + fnrt_c = currentCohort%prt%fnrt_c + store_c = currentCohort%prt%store_c ! Update the cohort's index within the size bin classes ! Update the cohort's index within the SCPF classification system @@ -902,7 +927,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentCohort => currentCohort%taller - enddo ! ends 'do while(associated(currentCohort)) + enddo co_loop1 ! ends 'do while(associated(currentCohort)) if ( currentPatch%total_canopy_area>currentPatch%area ) then if ( currentPatch%total_canopy_area-currentPatch%area > 0.001_r8 ) then @@ -915,15 +940,28 @@ subroutine canopy_summarization( nsites, sites, bc_in ) currentPatch%total_canopy_area = currentPatch%area endif + do ift=1,currentPatch%nupft + ft = currentPatch%unique_pfts(ift) + call UpdateSlowBiophysicalRates(currentPatch%bprates,ft,& + currentPatch%tveg_lpa%GetMean(), & + currentPatch%tveg_longterm%GetMean()) + end do + currentPatch => currentPatch%younger end do !patch loop - + call leaf_area_profile(sites(s)) if(hlm_radiation_model.eq.twostr_solver) then call FatesConstructRadElements(sites(s)) end if + currentPatch => sites(s)%oldest_patch + do while(associated(currentPatch)) + ! Transfer cohort info to patch level scratch space + call CopyCohortToCoArray(currentPatch) + currentPatch => currentPatch%younger + end do end do ! site loop return @@ -1037,6 +1075,8 @@ subroutine leaf_area_profile( currentSite ) ! or resized call cpatch%ReAllocateDynamics() + + ! These calls NaN and zero the above mentioned arrays call cpatch%NanDynamics() call cpatch%ZeroDynamics() @@ -1205,7 +1245,7 @@ subroutine leaf_area_profile( currentSite ) ! should have a value of exactly 1.0 in its top leaf layer ! -------------------------------------------------------------------------- - if ( (cpatch%NCL_p > 1) .and. & + if ( (cpatch%ncl > 1) .and. & (sum(cpatch%canopy_area_profile(1,:,1)) < 0.9999 )) then write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top' write(fates_log(), *) 'cl: ',1 @@ -1238,7 +1278,7 @@ subroutine leaf_area_profile( currentSite ) ! It should never be larger than 1 or less than 0. ! -------------------------------------------------------------------------- - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do iv = 1,cpatch%nleaf(cl,ft) if( debug .and. sum(cpatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then @@ -1298,7 +1338,7 @@ subroutine leaf_area_profile( currentSite ) ! preserve_b4b will be removed soon. This is kept here to prevent ! round off errors in the baseline tests for the two-stream code (RGK 12-27-23) if(preserve_b4b) then - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft do iv = 1, cpatch%nrad(cl,ft) if(cpatch%canopy_area_profile(cl,ft,iv) > 0._r8)then @@ -1308,16 +1348,24 @@ subroutine leaf_area_profile( currentSite ) end do end do else - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft if(cpatch%canopy_area_profile(cl,ft,1) > 0._r8 ) cpatch%canopy_mask(cl,ft) = 1 end do end do end if - end if if_any_canopy_area + ! Identify the maximum number of leaf layers for any canopy + ! or pft, this helps with efficient memory allocation during + ! photosynthesis + + cpatch%nleafmax = maxval(cpatch%nleaf) + + + + cpatch => cpatch%younger enddo !patch @@ -1574,14 +1622,14 @@ function calc_areaindex(cpatch,ai_type) result(ai) ai = 0._r8 if (trim(ai_type) == 'elai') then - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tlai') then - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tlai_profile(cl,ft,1:cpatch%nrad(cl,ft))) @@ -1589,14 +1637,14 @@ function calc_areaindex(cpatch,ai_type) result(ai) enddo elseif (trim(ai_type) == 'esai') then - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%esai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tsai') then - do cl = 1,cpatch%NCL_p + do cl = 1,cpatch%ncl do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tsai_profile(cl,ft,1:cpatch%nrad(cl,ft))) @@ -1711,7 +1759,7 @@ subroutine UpdateCohortLAI(currentCohort, canopy_layer_tlai, total_canopy_area) real(r8) :: treesai ! stem area index within crown m2/m2 ! Obtain the leaf carbon - leaf_c = currentCohort%prt%GetState(leaf_organ,carbon12_element) + leaf_c = sum(currentCohort%prt%leaf_c(:)) ! Note that tree_lai has an internal check on the canopy location call tree_lai_sai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & @@ -1754,4 +1802,188 @@ function NumCanopyLayers(currentPatch) result(z) end function NumCanopyLayers + ! ===================================================================================== + + subroutine CopyCohortToCoArray(patch) + + ! Copy data from the cohort linked-list to the cohort + ! vectors. These are only the variables that are + ! used at high frequencies + + type(fates_patch_type),intent(inout) :: patch + type(fates_cohort_type),pointer :: cohort + integer :: ico ! cohort index + integer :: ft + real(r8) :: leaf_c,leaf_n + real(r8) :: lnc_top + real(r8) :: store_c_target + real(r8) :: mr_reduction_factor + real(r8) :: storage_target_frac + real(r8) :: agb_frac ! fraction of biomass aboveground + real(r8) :: branch_frac ! fraction of aboveground woody biomass in branches + real(r8) :: crown_reduction ! reduction in crown biomass from damage + real(r8) :: sapw_c_bgw ! belowground sapwood + real(r8) :: sapw_c_agw ! aboveground sapwood + real(r8) :: sapw_c_undamaged ! the target sapwood of an undamaged tree + real(r8) :: sapw_n ! sapwood nitrogen + real(r8) :: sapw_n_bgw ! nitrogen in belowground portion of sapwood + real(r8) :: sapw_n_agw ! nitrogen in aboveground portion of sapwood + real(r8) :: sapw_n_undamaged ! nitrogen in sapwood of undamaged tree + + associate(coarr => patch%coarrays) + + ico = 0 + cohort => patch%tallest + do while (associated(cohort)) + ico = ico + 1 + + coarr%vcmax25top(ico) = cohort%vcmax25top + coarr%jmax25top(ico) = cohort%jmax25top + coarr%kp25top(ico) = cohort%kp25top + coarr%pft(ico) = cohort%pft + coarr%c_area(ico) = cohort%c_area + coarr%nplant(ico) = cohort%n + coarr%canopy_layer(ico) = cohort%canopy_layer + coarr%nv(ico) = cohort%nv + coarr%treesai(ico) = cohort%treesai + coarr%treelai(ico) = cohort%treelai + coarr%height(ico) = cohort%height + + call bleaf(cohort%dbh,cohort%pft,& + cohort%crowndamage,cohort%canopy_trim,1.0_r8,store_c_target) + + call storage_fraction_of_target(store_c_target, & + cohort%prt%store_c, & + storage_target_frac) + + call LowstorageMainRespReduction(storage_target_frac,cohort%pft, & + mr_reduction_factor) + + coarr%mr_reduction_factor(ico) = mr_reduction_factor + coarr%twostr_col(ico) = cohort%twostr_col + + ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) + ft = cohort%pft + select case(hlm_parteh_mode) + case (prt_carbon_allom_hyp) + lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/prt_params%slatop(ft) + case (prt_cnp_flex_allom_hyp) + leaf_c = sum(cohort%prt%leaf_c(:)) + if( (leaf_c*prt_params%slatop(ft)) > nearzero) then + leaf_n = sum(cohort%prt%leaf_n(:)) + lnc_top = leaf_n / (prt_params%slatop(ft) * leaf_c ) + else + lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/prt_params%slatop(ft) + end if + end select + + coarr%lnc_top(ico) = lnc_top + coarr%kn_leafn(ico) = DecayCoeffVcmax(cohort%vcmax25top, & + prt_params%leafn_vert_scaler_coeff1(ft), & + prt_params%leafn_vert_scaler_coeff2(ft)) + + coarr%kn_atk(ico) = DecayCoeffVcmax(cohort%vcmax25top, & + EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), & + EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft)) + + if (hlm_use_planthydro.eq.itrue ) then + coarr%leaf_psi(ico) = cohort%co_hydr%psi_ag(1) + coarr%btran(ico) = cohort%co_hydr%btran + else + coarr%btran(ico) = patch%btran_ft(ft) + coarr%leaf_psi(ico) = fates_unset_r8 + end if + + ! Organ Respiration: + ! -------------------------------------------------------------- + + if (hlm_use_tree_damage .eq. itrue) then + + ! Crown damage currenly only reduces the aboveground portion of + ! sapwood. Therefore we calculate the aboveground and the belowground portion + ! sapwood for use in stem respiration. + call GetCrownReduction(cohort%crowndamage, crown_reduction) + + else + crown_reduction = 0.0_r8 + end if + + ! If crown reduction is zero, undamaged sapwood target will equal sapwood carbon + agb_frac = prt_params%allom_agb_frac(cohort%pft) + branch_frac = param_derived%branch_frac(cohort%pft) + sapw_c_undamaged = cohort%prt%sapw_c / (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) + + ! Undamaged below ground portion + sapw_c_bgw = sapw_c_undamaged * (1.0_r8 - agb_frac) + + ! Damaged aboveground portion + sapw_c_agw = cohort%prt%sapw_c - sapw_c_bgw + + select case(hlm_parteh_mode) + case (prt_carbon_allom_hyp) + + coarr%live_stem_n(ico) = sapw_c_agw * & + prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + + coarr%live_croot_n(ico) = sapw_c_bgw * & + prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) + + coarr%fnrt_n(ico) = cohort%prt%fnrt_c * & + prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) + + case(prt_cnp_flex_allom_hyp) + + coarr%fnrt_n(ico) = cohort%prt%fnrt_n + + if (hlm_use_tree_damage .eq. itrue) then + + sapw_n_undamaged = cohort%prt%sapw_n / & + (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) + + sapw_n_bgw = sapw_n_undamaged * (1.0_r8 - agb_frac) + sapw_n_agw = cohort%prt%sapw_n - sapw_n_bgw + + coarr%live_croot_n(ico) = sapw_n_bgw + + coarr%live_stem_n(ico) = sapw_n_agw + else + + coarr%live_stem_n(ico) = prt_params%allom_agb_frac(cohort%pft) * & + cohort%prt%sapw_n + + coarr%live_croot_n(ico) = (1.0_r8-prt_params%allom_agb_frac(cohort%pft)) * & + cohort%prt%sapw_n + + end if + + case default + end select + + ! Force stem respiration to zero for non-woody species + if ( int(prt_params%woody(ft)) == ifalse) then + coarr%live_stem_n(ico) = 0._r8 + end if + + + ! IN/OUT (THIS SHOULDN'T BE NEEDED BUT BTRAN CALCULATION IS + ! OUT OF ORDER, THIS IS ONLY USED ONCE ON RESTARTS + coarr%g_sb_laweight(ico) = cohort%g_sb_laweight + + ! OUTPUTS + !coarr%resp_m_tstep(ico) = + !coarr%gpp_tstep(ico) = + !coarr%rdark_tstep(ico) = + !coarr%c13disc_clm(ico) = + !coarr%g_sb_laweight(ico) = + + + cohort => cohort%shorter + enddo + + coarr%ncohorts = ico + + end associate + + end subroutine CopyCohortToCoArray + end module EDCanopyStructureMod diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 016102dde0..f263abaf7b 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -2549,7 +2549,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in) height = EDPftvarcon_inst%hgt_min(ft) stem_drop_fraction = prt_params%phen_stem_drop_fraction(ft) fnrt_drop_fraction = prt_params%phen_fnrt_drop_fraction(ft) - l2fr = currentSite%rec_l2fr(ft, currentPatch%NCL_p) + l2fr = currentSite%rec_l2fr(ft, currentPatch%ncl) crowndamage = 1 ! new recruits are undamaged ! calculate DBH from initial height @@ -2782,7 +2782,7 @@ subroutine recruitment(currentSite, currentPatch, bc_in) call create_cohort(currentSite, currentPatch, ft, cohort_n, & height, 0.0_r8, dbh, prt, efleaf_coh, effnrt_coh, efstem_coh, & leaf_status, recruitstatus, init_recruit_trim, 0.0_r8, & - currentPatch%NCL_p, crowndamage, currentSite%spread, bc_in) + currentPatch%ncl, crowndamage, currentSite%spread, bc_in) ! Note that if hydraulics is on, the number of cohorts may have ! changed due to hydraulic constraints. @@ -3396,7 +3396,7 @@ subroutine UpdateRecruitStoich(csite) cpatch => csite%youngest_patch do while(associated(cpatch)) - cl = cpatch%ncl_p + cl = cpatch%ncl do ft = 1,numpft rec_l2fr_pft = csite%rec_l2fr(ft,cl) diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90 index 15195e4d0f..0a1501ed47 100644 --- a/biogeochem/FatesCohortMod.F90 +++ b/biogeochem/FatesCohortMod.F90 @@ -84,13 +84,12 @@ module FatesCohortMod real(r8) :: dbh ! diameter at breast height [cm] real(r8) :: coage ! age [years] real(r8) :: height ! height [m] - integer :: indexnumber ! unique number for each cohort (within clump?) integer :: canopy_layer ! canopy status of cohort [1 = canopy, 2 = understorey, etc.] real(r8) :: canopy_layer_yesterday ! recent canopy status of cohort [1 = canopy, 2 = understorey, etc.] ! real to be conservative during fusion integer :: crowndamage ! crown damage class of the cohort [1 = undamaged, >1 = damaged] real(r8) :: g_sb_laweight ! total conductance (stomata + boundary layer) of the cohort - ! weighted by its leaf area [m/s]*[m2] + ! weighted by its leaf area [m/s]*[m2] (NEEDED FOR RESTART) real(r8) :: canopy_trim ! fraction of the maximum leaf biomass that we are targeting [0-1] real(r8) :: leaf_cost ! how much does it cost to maintain leaves [kgC/m2/year] real(r8) :: excl_weight ! how much of this cohort is demoted each year, as a proportion of all cohorts @@ -141,19 +140,16 @@ module FatesCohortMod ! Units converted to a useful rate [kgC/indiv/yr] ! -------------------------------------------------------------------------- - real(r8) :: gpp_tstep ! Gross Primary Production (see above *) real(r8) :: gpp_acc real(r8) :: gpp_acc_hold real(r8) :: npp_acc real(r8) :: npp_acc_hold - real(r8) :: resp_m_tstep ! Maintenance respiration (see above *) real(r8) :: resp_m_acc real(r8) :: resp_m_acc_hold real(r8) :: resp_g_acc_hold - real(r8) :: c13disc_clm ! carbon 13 discrimination in new synthesized carbon at each indiv/timestep [ppm] real(r8) :: c13disc_acc ! carbon 13 discrimination in new synthesized carbon at each indiv/day ! at the end of a day [ppm] @@ -163,10 +159,8 @@ module FatesCohortMod real(r8) :: vcmax25top ! maximum carboxylation at canopy top and 25degC [umol CO2/m2/s] real(r8) :: jmax25top ! maximum electron transport rate at canopy top and 25degC [umol electrons/m2/s] - real(r8) :: tpu25top ! triose phosphate utilization rate at canopy top and 25degC [umol CO2/m2/s] real(r8) :: kp25top ! initial slope of CO2 response curve (C4 plants) at 25C - real(r8) :: ts_net_uptake(nlevleaf) ! net uptake of leaf layers [kgC/m2/timestep] real(r8) :: year_net_uptake(nlevleaf) ! net uptake of leaf layers [kgC/m2/year] ! used for CNP @@ -185,7 +179,6 @@ module FatesCohortMod ! in soil [kgN/plant/day] real(r8) :: sym_nfix_daily ! accumulated symbiotic N fixation from the roots [kgN/indiv/day] - real(r8) :: sym_nfix_tstep ! symbiotic N fixation from the roots for the time-step [kgN/indiv/timestep] real(r8) :: daily_n_gain ! sum of fixation and uptake of mineralized NH4/NO3 in solution as well as ! symbiotic fixation @@ -213,12 +206,6 @@ module FatesCohortMod ! to aid in reporting a more accurate sub-daily ! NEP - real(r8) :: rdark ! dark respiration [kgC/indiv/s] - real(r8) :: resp_m_unreduced ! diagnostic-only unreduced maintenance respiration [kgC/indiv/timestep] - real(r8) :: livestem_mr ! aboveground live stem maintenance respiration [kgC/indiv/s] - real(r8) :: livecroot_mr ! belowground live stem maintenance respiration [kgC/indiv/s] - real(r8) :: froot_mr ! live fine root maintenance respiration [kgC/indiv/s] - !--------------------------------------------------------------------------- ! DAMAGE @@ -361,7 +348,6 @@ subroutine NanValues(this) this%dbh = nan this%coage = nan this%height = nan - this%indexnumber = fates_unset_int this%canopy_layer = fates_unset_int this%canopy_layer_yesterday = nan this%crowndamage = fates_unset_int @@ -386,23 +372,18 @@ subroutine NanValues(this) this%size_class_lasttimestep = fates_unset_int ! CARBON AND NUTRIENT FLUXES - this%gpp_tstep = nan this%gpp_acc = nan this%gpp_acc_hold = nan this%npp_acc = nan this%npp_acc_hold = nan - this%resp_m_tstep = nan this%resp_m_acc = nan this%resp_m_acc_hold = nan this%resp_g_acc_hold = nan - this%c13disc_clm = nan this%c13disc_acc = nan this%vcmax25top = nan this%jmax25top = nan - this%tpu25top = nan this%kp25top = nan this%year_net_uptake(:) = nan - this%ts_net_uptake(:) = nan this%cnp_limiter = fates_unset_int this%cx_int = nan this%ema_dcxdt = nan @@ -412,7 +393,6 @@ subroutine NanValues(this) this%daily_nh4_uptake = nan this%daily_no3_uptake = nan this%sym_nfix_daily = nan - this%sym_nfix_tstep = nan this%daily_n_gain = nan this%daily_p_gain = nan this%daily_c_efflux = nan @@ -423,12 +403,7 @@ subroutine NanValues(this) this%seed_prod = nan ! RESPIRATION COMPONENTS - this%rdark = nan - this%resp_m_unreduced = nan this%resp_excess_hold = nan - this%livestem_mr = nan - this%livecroot_mr = nan - this%froot_mr = nan ! DAMAGE this%branch_frac = nan @@ -476,8 +451,6 @@ subroutine ZeroValues(this) ! ARGUMENTS class(fates_cohort_type), intent(inout) :: this - this%g_sb_laweight = 0._r8 - this%leaf_cost = 0._r8 this%excl_weight = 0._r8 this%prom_weight = 0._r8 @@ -490,12 +463,11 @@ subroutine ZeroValues(this) this%treesai = 0._r8 this%size_class = 1 this%coage_class = 1 - + + this%g_sb_laweight = 0._r8 this%size_class_lasttimestep = 0 - this%gpp_tstep = 0._r8 this%gpp_acc = 0._r8 this%npp_acc = 0._r8 - this%resp_m_tstep = 0._r8 this%resp_m_acc = 0._r8 ! do not zero these, they are not built @@ -507,10 +479,8 @@ subroutine ZeroValues(this) ! this%resp_g_acc_hold = nan ! this%resp_excess_hold = nan - this%c13disc_clm = 0._r8 this%c13disc_acc = 0._r8 - this%ts_net_uptake(:) = 0._r8 this%year_net_uptake(:) = 999._r8 ! this needs to be 999, or trimming of new cohorts will break. this%daily_nh4_uptake = 0._r8 @@ -533,12 +503,7 @@ subroutine ZeroValues(this) this%daily_n_demand = -9._r8 this%daily_p_demand = -9._r8 this%seed_prod = 0._r8 - this%rdark = 0._r8 - this%resp_m_unreduced = 0._r8 - this%livestem_mr = 0._r8 - this%livecroot_mr = 0._r8 - this%froot_mr = 0._r8 - + this%dmort = 0._r8 this%lmort_direct = 0._r8 this%lmort_collateral = 0._r8 @@ -684,8 +649,6 @@ subroutine Copy(this, copyCohort) class(fates_cohort_type), intent(in) :: this ! old cohort class(fates_cohort_type), intent(inout) :: copyCohort ! new cohort - copyCohort%indexnumber = fates_unset_int - ! POINTERS copyCohort%taller => NULL() copyCohort%shorter => NULL() @@ -724,22 +687,17 @@ subroutine Copy(this, copyCohort) copyCohort%size_class_lasttimestep = this%size_class_lasttimestep ! CARBON AND NUTRIENT FLUXES - copyCohort%gpp_tstep = this%gpp_tstep copyCohort%gpp_acc = this%gpp_acc copyCohort%gpp_acc_hold = this%gpp_acc_hold copyCohort%npp_acc = this%npp_acc copyCohort%npp_acc_hold = this%npp_acc_hold - copyCohort%resp_m_tstep = this%resp_m_tstep copyCohort%resp_m_acc = this%resp_m_acc copyCohort%resp_m_acc_hold = this%resp_m_acc_hold copyCohort%resp_g_acc_hold = this%resp_g_acc_hold - copyCohort%c13disc_clm = this%c13disc_clm copyCohort%c13disc_acc = this%c13disc_acc copyCohort%vcmax25top = this%vcmax25top copyCohort%jmax25top = this%jmax25top - copyCohort%tpu25top = this%tpu25top copyCohort%kp25top = this%kp25top - copyCohort%ts_net_uptake = this%ts_net_uptake copyCohort%year_net_uptake = this%year_net_uptake copyCohort%cnp_limiter = this%cnp_limiter @@ -753,7 +711,6 @@ subroutine Copy(this, copyCohort) copyCohort%daily_nh4_uptake = this%daily_nh4_uptake copyCohort%daily_no3_uptake = this%daily_no3_uptake copyCohort%sym_nfix_daily = this%sym_nfix_daily - copyCohort%sym_nfix_tstep = this%sym_nfix_tstep copyCohort%daily_n_gain = this%daily_n_gain copyCohort%daily_p_gain = this%daily_p_gain copyCohort%daily_c_efflux = this%daily_c_efflux @@ -764,12 +721,7 @@ subroutine Copy(this, copyCohort) copyCohort%seed_prod = this%seed_prod ! RESPIRATION COMPONENTS - copyCohort%rdark = this%rdark - copyCohort%resp_m_unreduced = this%resp_m_unreduced copyCohort%resp_excess_hold = this%resp_excess_hold - copyCohort%livestem_mr = this%livestem_mr - copyCohort%livecroot_mr = this%livecroot_mr - copyCohort%froot_mr = this%froot_mr ! DAMAGE copyCohort%branch_frac = this%branch_frac @@ -974,9 +926,6 @@ subroutine UpdateCohortBioPhysRates(this) this%jmax25top = sum(param_derived%jmax25top(ipft, 1:nleafage)* & frac_leaf_aclass(1:nleafage)) - this%tpu25top = sum(param_derived%tpu25top(ipft, 1:nleafage)* & - frac_leaf_aclass(1:nleafage)) - this%kp25top = sum(param_derived%kp25top(ipft, 1:nleafage)* & frac_leaf_aclass(1:nleafage)) @@ -984,14 +933,12 @@ subroutine UpdateCohortBioPhysRates(this) this%vcmax25top = EDPftvarcon_inst%vcmax25top(ipft, 1) this%jmax25top = param_derived%jmax25top(ipft, 1) - this%tpu25top = param_derived%tpu25top(ipft, 1) this%kp25top = param_derived%kp25top(ipft, 1) else this%vcmax25top = 0._r8 this%jmax25top = 0._r8 - this%tpu25top = 0._r8 this%kp25top = 0._r8 end if @@ -1088,7 +1035,7 @@ subroutine Dump(this) write(fates_log(),*) 'structural (dead) carbon = ', this%prt%GetState(struct_organ,carbon12_element) write(fates_log(),*) 'storage carbon = ', this%prt%GetState(store_organ,carbon12_element) write(fates_log(),*) 'reproductive carbon = ', this%prt%GetState(repro_organ,carbon12_element) - write(fates_log(),*) 'cohort%g_sb_laweight = ', this%g_sb_laweight + !write(fates_log(),*) 'cohort%g_sb_laweight = ', this%g_sb_laweight write(fates_log(),*) 'cohort%leaf_cost = ', this%leaf_cost write(fates_log(),*) 'cohort%canopy_layer = ', this%canopy_layer write(fates_log(),*) 'cohort%canopy_layer_yesterday = ', this%canopy_layer_yesterday @@ -1107,17 +1054,11 @@ subroutine Dump(this) write(fates_log(),*) 'cohort%coage_by_pft_class = ', this%coage_by_pft_class write(fates_log(),*) 'cohort%gpp_acc_hold = ', this%gpp_acc_hold write(fates_log(),*) 'cohort%gpp_acc = ', this%gpp_acc - write(fates_log(),*) 'cohort%gpp_tstep = ', this%gpp_tstep write(fates_log(),*) 'cohort%npp_acc_hold = ', this%npp_acc_hold write(fates_log(),*) 'cohort%npp_acc = ', this%npp_acc - write(fates_log(),*) 'cohort%resp_m_tstep = ', this%resp_m_tstep write(fates_log(),*) 'cohort%resp_m_acc = ', this%resp_m_acc write(fates_log(),*) 'cohort%resp_m_acc_hold = ', this%resp_m_acc_hold write(fates_log(),*) 'cohort%resp_g_acc_hold = ', this%resp_g_acc_hold - write(fates_log(),*) 'cohort%rdark = ', this%rdark - write(fates_log(),*) 'cohort%livestem_mr = ', this%livestem_mr - write(fates_log(),*) 'cohort%livecroot_mr = ', this%livecroot_mr - write(fates_log(),*) 'cohort%froot_mr = ', this%froot_mr write(fates_log(),*) 'cohort%dgmort = ', this%dgmort write(fates_log(),*) 'cohort%treelai = ', this%treelai write(fates_log(),*) 'cohort%treesai = ', this%treesai diff --git a/biogeochem/FatesPatchMod.F90 b/biogeochem/FatesPatchMod.F90 index 2f358cda63..096ca30d0a 100644 --- a/biogeochem/FatesPatchMod.F90 +++ b/biogeochem/FatesPatchMod.F90 @@ -34,7 +34,8 @@ module FatesPatchMod use FatesInterfaceTypesMod, only : numpft use shr_infnan_mod, only : nan => shr_infnan_nan, assignment(=) use shr_log_mod, only : errMsg => shr_log_errMsg - + use LeafBiophysicsMod, only : bprate_type + implicit none private @@ -60,6 +61,66 @@ module FatesPatchMod end type fates_cohort_vec_type + ! Scratch space for high frequency operations + ! in photosynthesis, respiration and + ! radiation scattering (soon) + + type :: fates_vecscr_type + + integer :: ncohorts ! Number of cohorts on this patch + + ! Inputs to photosynthesis + real(r8),allocatable :: vcmax25top(:) ! vcmax @ 25C top of canopy [umol/m2/s] + real(r8),allocatable :: jmax25top(:) ! jmax @ 25C top of canopy [umol/m2/s] + real(r8),allocatable :: kp25top(:) ! initial slope co2-curve @ 25C ... + integer,allocatable :: pft(:) + real(r8),allocatable :: c_area(:) ! areal extent of canopy [m2] + real(r8),allocatable :: nplant(:) ! number of plants per hectare + integer,allocatable :: canopy_layer(:) ! canopy index (1=top) + integer,allocatable :: nv(:) ! number of discrete vert veg layers + real(r8),allocatable :: treesai(:) ! stem area index of the plant [m2/m2] + real(r8),allocatable :: treelai(:) ! leaf area index of the plant [m2/m2] + real(r8),allocatable :: height(:) ! plant height [m] + real(r8),allocatable :: mr_reduction_factor(:) ! reduction factor for maint resp + real(r8),allocatable :: crown_reduction ! reduction in crown biomass from damage + real(r8),allocatable :: lnc_top(:) ! leaf nitrogen conc at the top [gn/m2] + real(r8),allocatable :: kn_leafn(:) ! leaf nitrogen vertical decay coeff + real(r8),allocatable :: kn_atk(:) ! atkin respiration vertical decay coeff + integer,allocatable :: twostr_col(:) ! this cohorts position index in the two-stream data + real(r8),allocatable :: btran(:) ! plant water stress response function [0-1] + + ! Inputs to organ respiration (fnrt,coarse root and above-ground stem) (i.e. nitrogen mass in kg) + real(r8),allocatable :: live_stem_n(ico) + real(r8),allocatable :: live_croot_n(ico) + real(r8),allocatable :: fnrt_n(ico) + + !Hydro only + real(r8),allocatable :: leaf_psi(:) !leaf suction [Mpa] + + ! Photosynthesis Outputs + real(r8),allocatable :: gpp_tstep(:) ! GPP [kgC/indiv/s] + real(r8),allocatable :: rdark_tstep(:) ! dark respiration [kgC/indiv/s] + real(r8),allocatable :: c13disc_clm(:) ! carbon 13 discrimination [ppm] + real(r8),allocatable :: g_sb_laweight(:) ! total conductance (stomata + + ! boundary layer) of the cohort + ! weighted by its leaf area [m/s]*[m2] + + ! Maintenance Respiration Outputs [kgC/indiv/s] + real(r8),allocatable :: livestem_mr(:) ! live stem MR + real(r8),allocatable :: froot_mr(:) ! fine root MR + real(r8),allocatable :: nfix_mr(:) ! respiration driving fixation + real(r8),allocatable :: sym_nfix_tstep(:) ! Symbiotic fixation rate kgN/plant/s + real(r8),allocatable :: livecroot_mr(:) ! live coarse root MR + real(r8),allocatable :: resp_m_tstep(:) ! total MR rate of plant (includes dark) + real(r8),allocatable :: resp_m_unreduced(:)! unreduced MR rate + + ! Cohort x leaf-layer + real(r8),allocatable :: ts_net_uptake(:,:) ! Net photosynthesis [kgC/m2/timestep] + + + + end type fates_vecscr_type + type, public :: fates_patch_type @@ -69,9 +130,15 @@ module FatesPatchMod type (fates_patch_type), pointer :: older => null() ! pointer to next older patch type (fates_patch_type), pointer :: younger => null() ! pointer to next younger patch type (fates_cohort_vec_type), pointer :: co_scr(:) ! Scratch vector of cohort properties - + + + ! Helper spaces for computational efficiency and pre-calcualtions outside + ! of high-frequency call spaces !--------------------------------------------------------------------------- + type(fates_vecscr_type),pointer :: coarrays ! Array structured cohort data + class(bprate_type),pointer :: bprates ! Biophysical rates + ! INDICES integer :: patchno ! unique number given to each new patch created for tracking integer :: nocomp_pft_label ! when nocomp is active, use this label for patch ID @@ -88,11 +155,19 @@ module FatesPatchMod integer :: age_class ! age class of the patch for history binning purposes real(r8) :: area ! patch area [m2] integer :: num_cohorts ! number of cohorts in patch - integer :: ncl_p ! number of occupied canopy layers + integer :: ncl ! number of occupied canopy layers + integer :: nleafmax ! maximum number of leaf layers required in + ! each of the canopy layers for evaluating photosynthesis + integer :: land_use_label ! patch label for land use classification (primaryland, secondaryland, etc) real(r8) :: age_since_anthro_disturbance ! average age for secondary forest since last anthropogenic disturbance [years] - logical :: changed_landuse_this_ts ! logical flag to track patches that have just undergone land use change [only used with nocomp and land use change] + logical :: changed_landuse_this_ts ! logical flag to track patches that have just undergone land use + ! change [only used with nocomp and land use change] + integer :: nupft ! number of unique PFTs in the patch + integer :: unique_pfts(maxpft) ! unique PFT indices that are on this patch + integer :: upft_index(maxpft) ! reverse of unique pfts + !--------------------------------------------------------------------------- ! RUNNING MEANS @@ -190,7 +265,7 @@ module FatesPatchMod ! ROOTS real(r8) :: btran_ft(maxpft) ! btran calculated seperately for each PFT real(r8) :: bstress_sal_ft(maxpft) ! bstress from salinity calculated seperately for each PFT - + !--------------------------------------------------------------------------- ! EXTERNAL SEED RAIN @@ -296,6 +371,13 @@ subroutine Init(this, num_swb, num_levsoil) allocate(this%sabs_dif(num_swb)) allocate(this%fragmentation_scaler(num_levsoil)) allocate(this%co_scr(max_cohort_per_patch)) + + ! Note that we allocate the innerds of + ! co-arrays during canopy summarization + allocate(this%coarrays) + + allocate(this%bprates) + call this%bprates%BPRateInitAlloc(numpft) ! initialize all values to nan call this%NanValues() @@ -329,8 +411,11 @@ subroutine ReAllocateDynamics(this) integer :: ncan ! Number of canopy layers integer :: prev_ncan ! Number of canopy layers previously ! as defined in the allocation space - - ncan = this%ncl_p + integer :: ncohorts + integer :: prev_ncohorts + type(fates_cohort_type),pointer :: cohort + + ncan = this%ncl nveg = maxval(this%nleaf(:,:)) ! Assume we will need to allocate, unless the @@ -382,26 +467,125 @@ subroutine ReAllocateDynamics(this) ! Add a little bit of buffer to the nveg ! so it doesn't need to be reallocated as much - nveg = nveg + 1 + + ! THESE SHOULD BE REARRANGED SO VEG IS INNER + ! THAT IS THE FAST ACCESS DIMENSION + allocate(this%tlai_profile(ncan,numpft,nveg+1)) + allocate(this%tsai_profile(ncan,numpft,nveg+1)) + allocate(this%elai_profile(ncan,numpft,nveg+1)) + allocate(this%esai_profile(ncan,numpft,nveg+1)) + allocate(this%canopy_area_profile(ncan,numpft,nveg+1)) + allocate(this%f_sun(ncan,numpft,nveg+1)) + allocate(this%fabd_sun_z(ncan,numpft,nveg+1)) + allocate(this%fabd_sha_z(ncan,numpft,nveg+1)) + allocate(this%fabi_sun_z(ncan,numpft,nveg+1)) + allocate(this%fabi_sha_z(ncan,numpft,nveg+1)) + allocate(this%nrmlzd_parprof_pft_dir_z(ncan,numpft,nveg+1)) + allocate(this%nrmlzd_parprof_pft_dif_z(ncan,numpft,nveg+1)) + allocate(this%ed_parsun_z(ncan,numpft,nveg+1)) + allocate(this%ed_parsha_z(ncan,numpft,nveg+1)) + allocate(this%ed_laisun_z(ncan,numpft,nveg+1)) + allocate(this%ed_laisha_z(ncan,numpft,nveg+1)) + end if + + ! Reallocate the Cohort Scratch array + ! --------------------------------------------------------------------------------- + + re_allocate = .true. + + ncohorts = 0 + cohort => this%tallest + do while (associated(cohort)) ! Cohort loop + ncohorts=ncohorts+1 + cohort => cohort%shorter + enddo - allocate(this%tlai_profile(ncan,numpft,nveg)) - allocate(this%tsai_profile(ncan,numpft,nveg)) - allocate(this%elai_profile(ncan,numpft,nveg)) - allocate(this%esai_profile(ncan,numpft,nveg)) - allocate(this%canopy_area_profile(ncan,numpft,nveg)) - allocate(this%f_sun(ncan,numpft,nveg)) - allocate(this%fabd_sun_z(ncan,numpft,nveg)) - allocate(this%fabd_sha_z(ncan,numpft,nveg)) - allocate(this%fabi_sun_z(ncan,numpft,nveg)) - allocate(this%fabi_sha_z(ncan,numpft,nveg)) - allocate(this%nrmlzd_parprof_pft_dir_z(ncan,numpft,nveg)) - allocate(this%nrmlzd_parprof_pft_dif_z(ncan,numpft,nveg)) - allocate(this%ed_parsun_z(ncan,numpft,nveg)) - allocate(this%ed_parsha_z(ncan,numpft,nveg)) - allocate(this%ed_laisun_z(ncan,numpft,nveg)) - allocate(this%ed_laisha_z(ncan,numpft,nveg)) + ! If the large patch arrays are not new, deallocate them + if(allocated(this%coarrays%vcmax25top)) then + prev_ncohorts = ubound(this%coarrays%vcmax25top,1) + prev_nveg = ubound(this%coarrays%ts_net_uptake,1) + if((ncohorts > prev_ncohorts) .or. (ncohorts < prev_ncohorts-3) .or. & + (nveg>prev_nveg) .or. (nveg ccohort%shorter + enddo + return + end subroutine FatesAccumulatePatchFluxes + + + end module EDAccumulateFluxesMod diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 177557aedf..aec9fb45a5 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -117,6 +117,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) integer :: j ! soil layer integer :: ifp ! patch vector index for the site integer :: ft ! plant functional type index + integer :: ico ! cohort index real(r8) :: smp_node ! matrix potential real(r8) :: rresis ! suction limitation to transpiration independent ! of root density @@ -193,9 +194,11 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) ! based on the cohort's leaf area. units: [m/s] * [m2] pftgs(1:maxpft) = 0._r8 + ico=0 ccohort => cpatch%tallest do while(associated(ccohort)) - pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%g_sb_laweight + ico=ico+1 + pftgs(ccohort%pft) = pftgs(ccohort%pft) + cpatch%coarrays%g_sb_laweight(ico) ccohort => ccohort%shorter enddo diff --git a/biogeophys/FatesLeafBiophysParamsMod.F90 b/biogeophys/FatesLeafBiophysParamsMod.F90 index fa578ad418..70b982f1b6 100644 --- a/biogeophys/FatesLeafBiophysParamsMod.F90 +++ b/biogeophys/FatesLeafBiophysParamsMod.F90 @@ -2,6 +2,8 @@ module FatesLeafBiophysParamsMod use FatesConstantsMod , only: r8 => fates_r8 use FatesConstantsMod , only: fates_check_param_set + use FatesConstantsMod , only: umolC_to_kgC + use FatesConstantsMod , only: g_per_kg use FatesGlobals, only : fates_log use FatesGlobals, only : endrun => fates_endrun use shr_log_mod , only : errMsg => shr_log_errMsg @@ -10,8 +12,6 @@ module FatesLeafBiophysParamsMod implicit none private ! Modules are private by default - save - public :: TransferParamsLeafBiophys public :: LeafBiophysReportParams @@ -58,6 +58,10 @@ subroutine TransferParamsLeafBiophys(pstruct) param_p => pstruct%GetParamFromName('fates_maintresp_leaf_ryan1991_baserate') allocate(lb_params%maintresp_leaf_ryan1991_baserate(numpft)) lb_params%maintresp_leaf_ryan1991_baserate(:) = param_p%r_data_1d(:) + + ! Convert from "gC/gN/s" to "umolC/gN/s @25 deg + lb_params%maintresp_leaf_ryan1991_baserate(:) = lb_params%maintresp_leaf_ryan1991_baserate(:) * & + (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) / (umolC_to_kgC * g_per_kg) param_p => pstruct%GetParamFromName('fates_maintresp_leaf_atkin2017_baserate') allocate(lb_params%maintresp_leaf_atkin2017_baserate(numpft)) @@ -102,10 +106,11 @@ subroutine TransferParamsLeafBiophys(pstruct) param_p => pstruct%GetParamFromName('fates_leaf_fnps') allocate(lb_params%fnps(numpft)) lb_params%fnps(:) = param_p%r_data_1d(:) - + return end subroutine TransferParamsLeafBiophys + ! ==================================================================================== subroutine LeafBiophysReportParams(is_master) @@ -117,8 +122,6 @@ subroutine LeafBiophysReportParams(is_master) character(len=32),parameter :: fmt_rout = '(a,F16.8)' character(len=32),parameter :: fmt_iout = '(a,I8)' - integer :: npft,ipft - if(debug_report .and. is_master) then write(fates_log(),fmt_iout) 'fates_leaf_c3psn = ',lb_params%c3psn write(fates_log(),fmt_iout) 'fates_leaf_stomatal_btran_model = ',lb_params%stomatal_btran_model diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90 index 1672ec5ff8..8221628fa8 100644 --- a/biogeophys/FatesPlantHydraulicsMod.F90 +++ b/biogeophys/FatesPlantHydraulicsMod.F90 @@ -2414,6 +2414,7 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) real(r8) :: sumcheck ! used to debug mass balance in soil horizon diagnostics integer :: nlevrhiz ! local for number of rhizosphere levels integer :: sc ! size class index + integer :: ico real(r8) :: lat,lon ! latitude and longitude of site real(r8) :: eff_por ! effective porosity real(r8) :: h2osoi_liqvol ! liquid water content [m3/m3] @@ -2508,11 +2509,13 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! ---------------------------------------------------------------------------- gscan_patch = 0.0_r8 + ico = 0 ccohort=>cpatch%tallest do while(associated(ccohort)) + ico=ico+1 ccohort_hydr => ccohort%co_hydr ccohort_hydr%psi_ag(1) = wrf_plant(leaf_p_media,ccohort%pft)%p%psi_from_th(ccohort_hydr%th_ag(1)) - gscan_patch = gscan_patch + ccohort%g_sb_laweight + gscan_patch = gscan_patch + cpatch%coarrays%g_sb_laweight(ico) ccohort => ccohort%shorter enddo !cohort @@ -2526,9 +2529,11 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + ico = 0 ccohort=>cpatch%tallest co_loop1: do while(associated(ccohort)) + ico = ico+1 ccohort_hydr => ccohort%co_hydr ft = ccohort%pft @@ -2537,9 +2542,9 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime) ! [mm H2O/plant/s] = [mm H2O/ m2 / s] * [m2 / patch] * [cohort/plant] * [patch/cohort] ! This can cause large transpiration due to small g_sb_laweight - if(ccohort%g_sb_laweight>nearzero) then + if(cpatch%coarrays%g_sb_laweight(ico)>nearzero) then qflx_tran_veg_indiv = bc_in(s)%qflx_transp_pa(ifp) * cpatch%total_canopy_area * & - (ccohort%g_sb_laweight/gscan_patch)/ccohort%n + (cpatch%coarrays%g_sb_laweight(ico)/gscan_patch)/ccohort%n else qflx_tran_veg_indiv = 0._r8 end if diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index 60cf234a64..192368a6dd 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -19,7 +19,24 @@ module FATESPlantRespPhotosynthMod ! ------------------------------------------------------------------------------------ ! !USES: - + use EDTypesMod , only : ed_site_type + use EDParamsMod , only : dinc_vai + use EDParamsMod , only : dlower_vai + use LeafBiophysicsMod , only : FastQ10 + use EDParamsMod , only : log_q10_mr_div10 + use FatesInterfaceTypesMod , only : bc_in_type + use FatesInterfaceTypesMod , only : bc_out_type + use EDCanopyStructureMod, only : calc_areaindex + use FatesConstantsMod, only : umolC_to_kgC + use FatesConstantsMod, only : umol_per_mmol + + use FatesParameterDerivedMod, only : param_derived + use FatesAllometryMod, only : bleaf, bstore_allom + use FatesAllometryMod, only : storage_fraction_of_target + use FatesAllometryMod, only : set_root_fraction + use DamageMainMod, only : GetCrownReduction + use FatesInterfaceTypesMod, only : hlm_use_tree_damage + use FatesGlobals, only : endrun => fates_endrun use FatesGlobals, only : fates_log use FatesGlobals, only : FatesWarn,N2S,A2S,I2S @@ -73,6 +90,11 @@ module FATESPlantRespPhotosynthMod use LeafBiophysicsMod, only : rsmax0 use LeafBiophysicsMod, only : DecayCoeffVcmax use LeafBiophysicsMod, only : VeloToMolarCF + use LeafBiophysicsMod, only : lb_params + use LeafBiophysicsMod, only : photosynth_acclim_model_none + use LeafBiophysicsMod, only : photosynth_acclim_model_kumarathunge_etal_2019 + use LeafBiophysicsMod, only : UpdateSlowBiophysicalRates + use LeafBiophysicsMod, only : UpdateFastBiophysicalRates use FatesRadiationMemMod, only : idirect ! CIME Globals @@ -82,45 +104,25 @@ module FATESPlantRespPhotosynthMod private public :: FatesPlantRespPhotosynthDrive ! Called by the HLM-Fates interface - + public :: FatesCondPhotoPatch + public :: FatesPlantRespPatch + character(len=*), parameter, private :: sourcefile = & __FILE__ - character(len=1024) :: warn_msg ! for defining a warning message logical :: debug = .false. + logical :: do_b4b = .false. !------------------------------------------------------------------------------------- contains !-------------------------------------------------------------------------------------- - subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) + subroutine FatesPlantRespPhotosynthDrive(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 EDTypesMod , only : ed_site_type - use EDParamsMod , only : dinc_vai - use EDParamsMod , only : dlower_vai - use FatesInterfaceTypesMod , only : bc_in_type - use FatesInterfaceTypesMod , only : bc_out_type - use EDCanopyStructureMod, only : calc_areaindex - use FatesConstantsMod, only : umolC_to_kgC - use FatesConstantsMod, only : umol_per_mmol - - use FatesParameterDerivedMod, only : param_derived - use FatesAllometryMod, only : bleaf, bstore_allom - use FatesAllometryMod, only : storage_fraction_of_target - use FatesAllometryMod, only : set_root_fraction - use DamageMainMod, only : GetCrownReduction - use FatesInterfaceTypesMod, only : hlm_use_tree_damage ! ARGUMENTS: ! ----------------------------------------------------------------------------------- @@ -130,14 +132,66 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) type(bc_out_type),intent(inout) :: bc_out(nsites) real(r8),intent(in) :: dtime + type (fates_patch_type) , pointer :: patch + integer :: s + + do s = 1,nsites + + patch => sites(s)%oldest_patch + do while (associated(patch)) + if(patch%nocomp_pft_label.ne.nocomp_bareground)then + + ! -------------------------------------------------------------------------- + ! 1 = initial value before the land-energy-balance solve + ! 2 = patch is currently marked for photosynthesis + ! 3 = patch has been called for photosynthesis already + ! -------------------------------------------------------------------------- + + if(bc_in(s)%filter_photo_pa(patch%patchno)==2)then + call FatesCondPhotoPatch(patch%patchno,sites(s),bc_in(s),bc_out(s)) + + call FatesPlantRespPatch(patch%patchno,sites(s),bc_in(s),dtime) + end if + end if + patch => patch%younger + end do + + end do + return + end subroutine FatesPlantRespPhotosynthDrive - ! LOCAL VARIABLES: + ! ================================================================================== + + subroutine FatesCondPhotoPatch(ifp,site,bc_in,bc_out) + + ! Refactor Objectives: + ! 1) Decrease memory footprint as much as possible + ! 2) Remove unecessary code + ! 3) Use better memory structure (remove heap memory) + ! 4) Evaluate loops to make sure no weird indexing + ! leaf layers are most accessed... + + ! Organizational notes: + ! 1) We do need to loop over cohorts, because we need to scale + ! up to patch level conductance. + ! 2) We could perform a lot of the filtering on the dyamics step + ! 3) These arrays for leaf area are not ideal... + + ! ----------------------------------------------------------------------------------- - type (fates_patch_type) , pointer :: currentPatch - type (fates_cohort_type), pointer :: currentCohort + ! !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 + ! ----------------------------------------------------------------------------------- + ! ----------------------------------------------------------------------------------- + integer,intent(in) :: ifp + type(ed_site_type),intent(inout) :: site + type(bc_in_type),intent(in) :: bc_in + type(bc_out_type),intent(inout) :: bc_out ! ----------------------------------------------------------------------------------- - ! These three arrays hold leaf-level biophysical rates that are calculated + ! The followingarrays 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 @@ -152,25 +206,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! are not modifying its order now. ! ----------------------------------------------------------------------------------- - ! leaf maintenance (dark) respiration [umol CO2/m**2/s] - real(r8) :: lmr_z(nlevleaf,maxpft,nclmax) + type (fates_patch_type), pointer :: patch + type (fates_cohort_type), pointer :: cohort - ! stomatal resistance [s/m] - real(r8) :: rs_z(nlevleaf,maxpft,nclmax) - - ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s] - real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax) - - ! photsynthesis - real(r8) :: psn_z(nlevleaf,maxpft,nclmax) - - ! carbon 13 in newly assimilated carbon at leaf level - real(r8) :: c13disc_z(nlevleaf,maxpft,nclmax) - ! Mask used to determine which leaf-layer biophysical rates have been - ! used already - logical :: rate_mask_z(nlevleaf,maxpft,nclmax) - 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 @@ -181,24 +220,12 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) 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) :: kn ! leaf nitrogen decay coefficient + real(r8) :: kn_leafn ! leaf nitrogen decay coefficient + real(r8) :: kn_atk ! maint resp decay coefficient for Atkin real(r8) :: gb_mol ! leaf boundary layer conductance (molar form: [umol /m**2/s]) real(r8) :: nscaler ! leaf nitrogen scaling coefficient real(r8) :: rdark_scaler ! scaling coefficient for Atkin dark respiration - real(r8) :: leaf_frac ! ratio of to leaf biomass to total alive biomass - real(r8) :: tcsoi ! Temperature response function for root respiration. - real(r8) :: tcwood ! Temperature response function for wood real(r8) :: patch_la ! exposed leaf area (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) :: sapw_c ! Sapwood carbon (kgC/plant) - real(r8) :: store_c_target ! Target storage carbon (kgC/plant) - real(r8) :: fnrt_c ! Fine root carbon (kgC/plant) - real(r8) :: fnrt_n ! Fine root nitrogen content (kgN/plant) - real(r8) :: leaf_c ! Leaf carbon (kgC/plant) - real(r8) :: leaf_n ! leaf nitrogen content (kgN/plant) real(r8) :: g_sb_leaves ! Mean combined (stomata+boundary layer) leaf conductance [m/s] ! over all of the patch's leaves. The "sb" refers to the combined ! "s"tomatal and "b"oundary layer. @@ -209,49 +236,23 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: r_stomata ! Mean stomatal resistance across all leaves in the patch [s/m] real(r8) :: maintresp_reduction_factor ! factor by which to reduce maintenance ! respiration when storage pools are low - real(r8) :: b_leaf ! leaf biomass kgC - real(r8) :: frac ! storage pool as a fraction of target leaf biomass - ! over each cohort x layer. + real(r8) :: storage_target_frac ! carbon storage as a fraciton of the target value + real(r8) :: store_c_target ! Target storage carbon (kgC/plant) real(r8) :: cohort_eleaf_area ! This is the effective leaf area [m2] reported by each cohort real(r8) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] - real(r8) :: lmr25top ! canopy top leaf maint resp rate at 25C - ! for this plant or pft (umol CO2/m**2/s) real(r8) :: lai_canopy_above ! the LAI in the canopy layers above the layer of interest real(r8) :: leaf_veg_frac ! fraction of vegetation area (leaf+stem) that is just leaf real(r8) :: cumulative_lai ! the cumulative LAI, top down, to the leaf layer of interest real(r8) :: leaf_psi ! leaf xylem matric potential [MPa] (only meaningful/used w/ hydro) - real(r8) :: fnrt_mr_layer ! fine root maintenance respiation per layer [kgC/plant/s] + integer :: isunsha ! Index for differentiating sun and shade - real(r8) :: fnrt_mr_nfix_layer ! fineroot maintenance respiration - ! specifically for symbiotic fixation [kgC/plant/layer/s] - real(r8) :: nfix_layer ! Nitrogen fixed in each layer this timestep [kgN/plant/layer/timestep] - real(r8), allocatable :: rootfr_ft(:,:) ! Root fractions per depth and PFT - real(r8) :: agb_frac ! fraction of biomass aboveground - real(r8) :: branch_frac ! fraction of aboveground woody biomass in branches - real(r8) :: crown_reduction ! reduction in crown biomass from damage - real(r8) :: sapw_c_bgw ! belowground sapwood - real(r8) :: sapw_c_agw ! aboveground sapwood - real(r8) :: sapw_c_undamaged ! the target sapwood of an undamaged tree - real(r8) :: sapw_n ! sapwood nitrogen - real(r8) :: sapw_n_bgw ! nitrogen in belowground portion of sapwood - real(r8) :: sapw_n_agw ! nitrogen in aboveground portion of sapwood - real(r8) :: sapw_n_undamaged ! nitrogen in sapwood of undamaged tree real(r8) :: rd_abs_leaf, rb_abs_leaf ! Watts of diffuse and beam light absorbed by leaves over this ! depth interval and ground footprint (m2) - real(r8) :: r_abs_stem, r_abs_snow ! Watts of light absorbed by stems and snow over this depth interval and - ! ground footprint real(r8) :: rb_abs, rd_abs ! Total beam and diffuse radiation absorbed over this depth interval ! and ground footprint real(r8) :: fsun ! sun-shade fraction real(r8) :: par_per_sunla, par_per_shala ! PAR per sunlit and shaded leaf area [W/m2 leaf] - real(r8) :: ac_utest, aj_utest ! Gross rubisco and rubp limited assimilation (for unit tests) - real(r8) :: ap_utest, co2_inter_c_utest ! PEP limited assimilation, and intracellular co2 (for unit tests) - real(r8),dimension(150) :: cohort_vaitop ! The top-down integrated vegetation area index - ! (leaf+stem) at the top of the layer - real(r8),dimension(150) :: cohort_vaibot ! The top-down integrated vegetation area index - ! (leaf+stem) at the bottom of the layer - real(r8),dimension(150) :: cohort_layer_elai ! exposed leaf area index of the layer - real(r8),dimension(150) :: cohort_layer_esai ! exposed stem area index of the layer + real(r8) :: co2_inter_c ! intracellular co2 (for unit tests) (Pa) real(r8) :: cohort_elai ! exposed leaf area index of the cohort real(r8) :: cohort_esai ! exposed stem area index of the cohort real(r8) :: laisun,laisha ! m2 of exposed or shaded leaves per m2 of crown @@ -260,7 +261,6 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: canopy_area ! For Norman radiation, the fraction of crown area per ! canopy area in each layer ! they will sum to 1.0 in the fully closed canopy layers - real(r8) :: elai ! effective leaf area index real(r8) :: area_frac ! this is either f_sun , or 1-f_sun, its for area weighting real(r8) :: psn_ll ! Leaf level photosyntheis real(r8) :: gstoma_ll ! leaf level stomatal conductance (separate for sun shade) [m/s] @@ -273,884 +273,769 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) real(r8) :: gs2 ! optional btran scaling factor for Medlyn conductance only, instead ! of applying to gs1, this would scale the whole non-intercept ! portion of the conductance equation - real(r8) :: vmol_cf ! velocity to molar conductance conversion (m/s) -> (umol/m2/s) + real(r8) :: leaf_c ! Leaf carbon (kgC/plant) + real(r8) :: leaf_n ! leaf nitrogen content (kgN/plant) + integer :: cl,iv,ft,uft,ift,ico ! indices + integer :: nv ! number of leaf layers + integer :: solve_iter ! number of iterations required for photosynthesis solve + + ! Temporary for comparing radiation solutions + real(r8) :: Rb_abs_el,Rd_abs_el,rd_abs_leaf_el,rb_abs_leaf_el,r_abs_stem_el,r_abs_snow_el,leaf_sun_frac + logical :: call_fail + - ! ----------------------------------------------------------------------------------- - ! 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 - integer :: NCL_p ! number of canopy layers in patch - integer :: iage ! loop counter for leaf age classes - integer :: solve_iter ! number of iterations required for photosynthesis solve - - ! Parameters + ! Constants ! Absolute convergence tolerance on solving intracellular CO2 concentration [Pa] real(r8), parameter :: ci_tol = 0.5_r8 - - ! Base rate is at 20C. Adjust to 25C using the CN Q10 = 1.5 - ! (gC/gN/s) - ! ------------------------------------------------------------------------ + patch => site%pa_vec(ifp)%p - ! ----------------------------------------------------------------------------------- - ! Photosynthesis and stomatal conductance parameters, from: - ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 - ! ----------------------------------------------------------------------------------- + if(.not.associated(patch))then + write(fates_log(),*)'did not find ifp?' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + bc_out%rssun_pa(ifp) = 0._r8 + bc_out%rssha_pa(ifp) = 0._r8 + patch%c_lblayer = fates_unset_r8 + if(patch%num_cohorts == 0) return + + associate(co_arr => patch%coarrays) + + ! We use this block to enable the use of automatic + ! arrays (stack) as opposed to dynamic allocations + ! (heap), which enables faster computation + block + + ! leaf maintenance (dark) respiration [umol CO2/m**2/s] + real(r8) :: lmr_z(patch%nleafmax,patch%nupft,patch%ncl) + + ! stomatal resistance [s/m] + real(r8) :: rs_z(patch%nleafmax,patch%nupft,patch%ncl) + + ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s] + real(r8) :: anet_av_z(patch%nleafmax,patch%nupft,patch%ncl) + + ! photsynthesis + real(r8) :: psn_z(patch%nleafmax,patch%nupft,patch%ncl) + + ! carbon 13 in newly assimilated carbon at leaf level + real(r8) :: c13disc_z(patch%nleafmax,patch%nupft,patch%ncl) + + ! Mask used to determine which leaf-layer + ! biophysical rates have been used already + logical :: rate_mask_z(patch%nleafmax,patch%nupft,patch%ncl) + + real(r8) :: cohort_vaitop(patch%nleafmax) ! The top-down integrated vegetation area index + ! (leaf+stem) at the top of the layer + real(r8) :: cohort_vaibot(patch%nleafmax) ! The top-down integrated vegetation area index + ! (leaf+stem) at the bottom of the layer + real(r8) :: cohort_layer_elai(patch%nleafmax) ! exposed leaf area index of the layer + real(r8) :: cohort_layer_esai(patch%nleafmax) ! exposed stem area index of the layer + + g_sb_leaves = 0._r8 + patch_la = 0._r8 + + ! Update temperature scaling functions on biophysical rates + ! we only need to do this once per PFT + do ift=1,patch%nupft + ft = patch%unique_pfts(ift) + if(do_b4b)then + call UpdateSlowBiophysicalRates(patch%bprates,ft, & + patch%tveg_lpa%GetMean(), & + patch%tveg_longterm%GetMean()) + end if + call UpdateFastBiophysicalRates(patch%bprates,ft, & + bc_in%t_veg_pa(ifp),lb_params%c3psn(ft),hlm_maintresp_leaf_model) + end do + + ! Part II. 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) + ! leaf boundary layer conductance of h20 + ! constrained vapor pressure + + call GetCanopyGasParameters(bc_in%forc_pbot, & ! in + bc_in%oair_pa(ifp), & ! in + bc_in%t_veg_pa(ifp), & ! in + mm_kco2, & ! out + mm_ko2, & ! out + co2_cpoint) + + ! The host models use velocity based conductances and resistance + ! this is the factor that converts a conductance from + ! [m/s] to [umol/m2/s] + vmol_cf = VeloToMolarCF(bc_in%forc_pbot,bc_in%tgcm_pa(ifp)) + + ! ------------------------------------------------------------------------ + ! 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. + ! NOTE: Only need to flush mask on the number of used pfts, not the whole + ! scratch space. + ! ------------------------------------------------------------------------ + rate_mask_z(:,:,:) = .false. + + do ico = 1,co_arr%ncohorts + + ! Identify the canopy layer (cl), functional type (ft) + ! and the leaf layer (IV) for this cohort + ft = co_arr%pft(ico) + uft = patch%upft_index(ft) !unique pft index (1:npfts_unique) + cl = co_arr%canopy_layer(ico) + + ! Calculate the cohort specific elai profile + ! And the top and bottom edges of the veg area index + ! of each layer bin are. Note, if the layers + ! sink below the ground snow line, then the effective + ! LAI and SAI start to shrink to zero, as well as + ! the difference between vaitop and vaibot. + if(co_arr%treesai(ico)>0._r8)then + do iv = 1,co_arr%nv(ico) + call VegAreaLayer(co_arr%treelai(ico), & + co_arr%treesai(ico), & + co_arr%height(ico), & + iv, & + co_arr%nv(ico), & + co_arr%pft(ico), & + site%snow_depth, & + cohort_vaitop(iv), & + cohort_vaibot(iv), & + cohort_layer_elai(iv), & + cohort_layer_esai(iv)) + end do + cohort_elai = sum(cohort_layer_elai(1:co_arr%nv(ico))) + cohort_esai = sum(cohort_layer_esai(1:co_arr%nv(ico))) + else + cohort_vaitop(:) = 0._r8 + cohort_vaibot(:) = 0._r8 + cohort_elai = 0._r8 + cohort_esai = 0._r8 + end if + + !co_arr%mr_reduction_factor(ico) + + ! are there any leaves of this pft in this layer? + canopy_mask_if: if(patch%canopy_mask(cl,ft) == 1)then + + ! Loop over leaf-layers + + leaf_veg_frac = co_arr%treelai(ico)/(co_arr%treelai(ico)+co_arr%treesai(ico)) + + leaf_layer_loop : do iv = 1,co_arr%nv(ico) + + ! ------------------------------------------------------------ + ! 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. + ! Other cases where we need to solve for every cohort + ! in every leaf layer: nutrient dynamic mode, multiple leaf + ! age classes + ! ------------------------------------------------------------ + + rate_mask_if: if ( .not.rate_mask_z(iv,uft,cl) .or. & + (hlm_use_planthydro.eq.itrue) .or. & + (hlm_radiation_model .eq. twostr_solver ) .or. & + (nleafage > 1) .or. & + (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then + + ! These values are incremented, therefore since + ! sometimes we re-do layers, we need to re-zero them as well + ! since it is not an over-write + psn_z(iv,uft,cl) = 0._r8 + anet_av_z(iv,uft,cl) = 0._r8 + c13disc_z(iv,uft,cl) = 0._r8 + + if (hlm_use_planthydro.eq.itrue ) then + + btran_eff = co_arr%btran(ico) + + ! Find the cumulative LAI from the top of this cohort's crown to the + ! center of the current veg layer. If this is the cohort's last layer + ! then the mid-point is between the dlower and the total lai + + lai_canopy_above = sum(patch%canopy_layer_tlai(1:cl-1)) + + if(iv == co_arr%nv(ico)) then + cumulative_lai = lai_canopy_above + leaf_veg_frac * (dlower_vai(iv)+0.5_r8*(co_arr%treelai(ico)+co_arr%treesai(ico)-dlower_vai(iv))) + else + cumulative_lai = lai_canopy_above + leaf_veg_frac * (dlower_vai(iv)+0.5_r8*dinc_vai(iv)) + end if + + leaf_psi = co_arr%leaf_psi(ico) + + else + + btran_eff = patch%btran_ft(ft) + ! For consistency sake, we use total LAI here, and not exposed + ! if the plant is under-snow, it will be effectively dormant for + ! the purposes of nscaler + + cumulative_lai = sum(patch%canopy_layer_tlai(1:cl-1)) + & + sum(patch%tlai_profile(cl,ft,1:iv-1)) + & + 0.5*patch%tlai_profile(cl,ft,iv) + + leaf_psi = fates_unset_r8 + + end if + + if(do_fates_salinity)then + btran_eff = btran_eff*patch%bstress_sal_ft(ft) + endif + + ! Scale for leaf nitrogen profile + nscaler = exp(-co_arr%kn_leafn(ico) * cumulative_lai) + + ! 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 + + ! Part VII: Calculate dark respiration (leaf maintenance) for this layer + + select case (hlm_maintresp_leaf_model) + + case (lmrmodel_ryan_1991) + + call LeafLayerMaintenanceRespiration_Ryan_1991( co_arr%lnc_top(ico), & ! in + nscaler, & ! in + ft, & ! in + bc_in%t_veg_pa(ifp), & ! in + patch%bprates, & ! in + lmr_z(iv,uft,cl)) ! out + case (lmrmodel_atkin_etal_2017) - associate( & - slatop => prt_params%slatop , & ! specific leaf area at top of canopy, - ! projected area basis [m^2/gC] - woody => prt_params%woody) ! Is vegetation woody or not? - - 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 - - - - ! Pre-process some variables that are PFT dependent - ! but not environmentally dependent - ! ------------------------------------------------------------------------ - - allocate(rootfr_ft(numpft, bc_in(s)%nlevsoil)) - - do ft = 1,numpft - call set_root_fraction(rootfr_ft(ft,:), ft, & - bc_in(s)%zi_sisl, & - bc_in(s)%max_rooting_depth_index_col) - end do - - - currentpatch => sites(s)%oldest_patch - do while (associated(currentpatch)) - - ifp = currentPatch%patchno - - if_notbare: if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then - - NCL_p = currentPatch%NCL_p - - ! Part I. Zero output boundary conditions - ! --------------------------------------------------------------------------- - bc_out(s)%rssun_pa(ifp) = 0._r8 - bc_out(s)%rssha_pa(ifp) = 0._r8 - - g_sb_leaves = 0._r8 - patch_la = 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 already - ! --------------------------------------------------------------------------- - if_filter2: 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: - ! currentPatch%ncan(:,:) - ! currentPatch%canopy_mask(:,:) - call UpdateCanopyNCanNRadPresent(currentPatch) - - - ! 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) - ! 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 - mm_kco2, & ! out - mm_ko2, & ! out - co2_cpoint) - - ! The host models use velocity based conductances and resistance - ! this is the factor that converts a conductance from - ! [m/s] to [umol/m2/s] - vmol_cf = VeloToMolarCF(bc_in(s)%forc_pbot,bc_in(s)%tgcm_pa(ifp)) - - ! ------------------------------------------------------------------------ - ! 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. - ! NOTE: Only need to flush mask on the number of used pfts, not the whole - ! scratch space. - ! ------------------------------------------------------------------------ - rate_mask_z(:,1:numpft,:) = .false. - psn_z(:,:,:) = 0._r8 - anet_av_z(:,:,:) = 0._r8 - c13disc_z(:,:,:) = 0._r8 - rs_z(:,:,:) = 0._r8 - lmr_z(:,:,:) = 0._r8 - - if_any_cohorts: if(currentPatch%num_cohorts > 0)then - - currentCohort => currentPatch%tallest - do_cohort_drive: 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 - - ! Calculate the cohort specific elai profile - ! And the top and bottom edges of the veg area index - ! of each layer bin are. Note, if the layers - ! sink below the ground snow line, then the effective - ! LAI and SAI start to shrink to zero, as well as - ! the difference between vaitop and vaibot. - if(currentCohort%treesai>0._r8)then - do iv = 1,currentCohort%nv - call VegAreaLayer(currentCohort%treelai, & - currentCohort%treesai, & - currentCohort%height, & - iv, & - currentCohort%nv, & - currentCohort%pft, & - sites(s)%snow_depth, & - cohort_vaitop(iv), & - cohort_vaibot(iv), & - cohort_layer_elai(iv), & - cohort_layer_esai(iv)) - end do - - cohort_elai = sum(cohort_layer_elai(1:currentCohort%nv)) - cohort_esai = sum(cohort_layer_esai(1:currentCohort%nv)) + ! This uses the relationship between leaf N and respiration from Atkin et al + ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler. + ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the + ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration + ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis + ! this number can be smaller. There is some observational evidence for this being the case + ! in Lamour et al. 2023. + + rdark_scaler = exp(-co_arr%kn_atk(ico) * cumulative_lai) + call LeafLayerMaintenanceRespiration_Atkin_etal_2017(co_arr%lnc_top(ico), & ! in + rdark_scaler, & ! in + ft, & ! in + bc_in%t_veg_pa(ifp), & ! in + patch%tveg_lpa%GetMean(), & ! in + patch%bprates, & ! in + lmr_z(iv,uft,cl)) ! out - else - cohort_layer_elai(:) = 0._r8 - cohort_layer_esai(:) = 0._r8 - cohort_vaitop(:) = 0._r8 - cohort_vaibot(:) = 0._r8 - cohort_elai = 0._r8 - cohort_esai = 0._r8 - end if + case default - ! MLO. Assuming target to be related to leaf biomass when leaves are fully - ! flushed. But unsure whether this call is correct or not, shouldn't we get - ! the target value directly from the bstore_allom? - call bleaf(currentCohort%dbh,currentCohort%pft,& - currentCohort%crowndamage,currentCohort%canopy_trim,1.0_r8,store_c_target) - ! call bstore_allom(currentCohort%dbh,currentCohort%pft, & - ! currentCohort%canopy_trim,store_c_target) - - call storage_fraction_of_target(store_c_target, & - currentCohort%prt%GetState(store_organ, carbon12_element), & - frac) - call LowstorageMainRespReduction(frac,currentCohort%pft, & - maintresp_reduction_factor) - - ! are there any leaves of this pft in this layer? - canopy_mask_if: if(currentPatch%canopy_mask(cl,ft) == 1)then - - ! Loop over leaf-layers - - leaf_veg_frac = currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai) - - leaf_layer_loop : 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), - ! 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. - ! Other cases where we need to solve for every cohort - ! in every leaf layer: nutrient dynamic mode, multiple leaf - ! age classes - ! ------------------------------------------------------------ - - rate_mask_if: if ( .not.rate_mask_z(iv,ft,cl) .or. & - (hlm_use_planthydro.eq.itrue) .or. & - (hlm_radiation_model .eq. twostr_solver ) .or. & - (nleafage > 1) .or. & - (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then - - - ! These values are incremented, therefore since - ! sometimes we re-do layers, we need to re-zero them as well - ! since it is not an over-write - psn_z(iv,ft,cl) = 0._r8 - anet_av_z(iv,ft,cl) = 0._r8 - c13disc_z(iv,ft,cl) = 0._r8 - - if (hlm_use_planthydro.eq.itrue ) then - - btran_eff = currentCohort%co_hydr%btran - - ! Find the cumulative LAI from the top of this cohort's crown to the - ! center of the current veg layer. If this is the cohort's last layer - ! then the mid-point is between the dlower and the total lai - - lai_canopy_above = sum(currentPatch%canopy_layer_tlai(1:cl-1)) - - if(iv == currentCohort%nv) then - cumulative_lai = lai_canopy_above + leaf_veg_frac * (dlower_vai(iv)+0.5_r8*(currentCohort%treelai+currentCohort%treesai-dlower_vai(iv))) - else - cumulative_lai = lai_canopy_above + leaf_veg_frac * (dlower_vai(iv)+0.5_r8*dinc_vai(iv)) - end if - - leaf_psi = currentCohort%co_hydr%psi_ag(1) - - else - - btran_eff = currentPatch%btran_ft(ft) - ! For consistency sake, we use total LAI here, and not exposed - ! if the plant is under-snow, it will be effectively dormant for - ! the purposes of nscaler - - cumulative_lai = sum(currentPatch%canopy_layer_tlai(1:cl-1)) + & - sum(currentPatch%tlai_profile(cl,ft,1:iv-1)) + & - 0.5*currentPatch%tlai_profile(cl,ft,iv) - - leaf_psi = fates_unset_r8 - - end if - - if(do_fates_salinity)then - btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft) - endif - - ! 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 - - kn = DecayCoeffVcmax(currentCohort%vcmax25top, & - prt_params%leafn_vert_scaler_coeff1(ft), & - prt_params%leafn_vert_scaler_coeff2(ft)) - - ! Scale for leaf nitrogen profile - nscaler = exp(-kn * cumulative_lai) - - ! 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 - ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf) - select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) - - lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - - case (prt_cnp_flex_allom_hyp) - - leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) - if( (leaf_c*slatop(ft)) > nearzero) then - leaf_n = currentCohort%prt%GetState(leaf_organ, nitrogen_element) - lnc_top = leaf_n / (slatop(ft) * leaf_c ) - else - lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - end if - - ! If one wants to break coupling with dynamic N conentrations, - ! use the stoichiometry parameter - ! lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) - - end select - - - ! Part VII: Calculate dark respiration (leaf maintenance) for this layer - - select case (hlm_maintresp_leaf_model) - - case (lmrmodel_ryan_1991) - - call LeafLayerMaintenanceRespiration_Ryan_1991( lnc_top, & ! in - nscaler, & ! in - ft, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - lmr_z(iv,ft,cl)) ! out - - case (lmrmodel_atkin_etal_2017) - - ! This uses the relationship between leaf N and respiration from Atkin et al - ! for the top of the canopy, but then scales through the canopy based on a rdark_scaler. - ! To assume proportionality with N through the canopy following Lloyd et al. 2010, use the - ! default parameter value of 2.43, which results in the scaling of photosynthesis and respiration - ! being proportional through the canopy. To have a steeper decrease in respiration than photosynthesis - ! this number can be smaller. There is some observational evidence for this being the case - ! in Lamour et al. 2023. - - kn = DecayCoeffVcmax(currentCohort%vcmax25top, & - EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1(ft), & - EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2(ft)) - - rdark_scaler = exp(-kn * cumulative_lai) - - call LeafLayerMaintenanceRespiration_Atkin_etal_2017( lnc_top, & ! in - rdark_scaler, & ! in - ft, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - currentPatch%tveg_lpa%GetMean(), & ! in - lmr_z(iv,ft,cl)) ! out - - case default - - write (fates_log(),*)'error, incorrect leaf respiration model specified' - call endrun(msg=errMsg(sourcefile, __LINE__)) - - end select - - ! Pre-process PAR absorbed per unit leaf area for different schemes - ! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves] - ! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves] - ! fsun = [m2 of sunlit leaves / m2 of total leaves] - ! laisun: m2 of exposed leaf, per m2 of crown. If this is the lowest layer - ! for the pft/canopy group, than the m2 per crown is probably not - ! as large as the layer above. - ! ------------------------------------------------------------------ - - if_radsolver: if(hlm_radiation_model.eq.norman_solver) then - - laisun = currentPatch%ed_laisun_z(cl,ft,iv) - laisha = currentPatch%ed_laisha_z(cl,ft,iv) - par_per_sunla = currentPatch%ed_parsun_z(cl,ft,iv) - par_per_shala = currentPatch%ed_parsha_z(cl,ft,iv) - canopy_area = currentPatch%canopy_area_profile(cl,ft,iv) - fsun = currentPatch%f_sun(cl,ft,iv) - - else ! Two-stream - - if(cohort_layer_elai(iv) > nearzero .and. sites(s)%coszen>0._r8 ) then - - call FatesGetCohortAbsRad(currentPatch, currentCohort, ipar, & - cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, & - rb_abs, rd_abs, rb_abs_leaf, rd_abs_leaf, fsun) - - ! rd_abs_leaf: Watts of diffuse light absorbed by leaves over this - ! depth interval and ground footprint (m2) - ! rd_abs_leaf*fsun Watts of diffuse light absorbed by sunlit leaves - ! over this depth interval and ground footprint (m2) - ! rb_abs_leaf Watts of beam absorbed by sunlit leaves over this - ! depth interval and ground footprint (m2) - ! cohort_layer_elai*fsun Leaf area in sunlight within this interval and ground footprint - ! cohort_layer_elai*(1-fsun) Leaf area in shade within this interval and ground footprint - - laisun = (fsun*cohort_layer_elai(iv)) - laisha = ((1._r8 - fsun)*cohort_layer_elai(iv)) - if(fsun>nearzero) then - par_per_sunla = (rd_abs_leaf*fsun + rb_abs_leaf)! / laisun - else - par_per_sunla = 0._r8 - end if - par_per_shala = rd_abs_leaf*(1._r8-fsun) !/ laisha - canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) - - else - - par_per_sunla = 0._r8 - par_per_shala = 0._r8 - laisun = 0.5_r8*cohort_layer_elai(iv) - laisha = 0.5_r8*cohort_layer_elai(iv) - canopy_area = 1._r8 !currentPatch%canopy_area_profile(cl,ft,iv) - fsun = 0.5_r8 !avoid div0, should have no impact - - end if - - end if if_radsolver - - ! Perform photosynthesis calculations on sunlit and shaded leaves - ! --------------------------------------------------------------- - - ! Calculate leaf boundary layer conductance in molar form [umol/m2/s] - gb_mol = (1._r8/bc_in(s)%rb_pa(ifp)) * vmol_cf - - gstoma = 0._r8 - do_sunsha: do isunsha = 1,2 - - ! Determine absorbed PAR per square meter of leaf - ! If there is no leaf area perform a trivial solution - - if(isunsha == idirect) then - leaf_area = laisun*canopy_area - par_abs = ConvertPar(leaf_area, par_per_sunla) - area_frac = fsun - else - leaf_area = laisha*canopy_area - par_abs = ConvertPar(leaf_area, par_per_shala) - area_frac = 1._r8 - fsun - end if - - if( leaf_area < nearzero ) then - - ! Note: With no leaf area do not increment - ! any fluxes. Assume a nominal conductance - ! of maximum resistance - gstoma = gstoma + area_frac/rsmax0 - - cycle do_sunsha - end if - - ! 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(ft, & ! in - currentCohort%vcmax25top, & ! in - currentCohort%jmax25top, & ! in - currentCohort%kp25top, & ! in - nscaler, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%dayl_factor_pa(ifp), & ! in - currentPatch%tveg_lpa%GetMean(), & ! in - currentPatch%tveg_longterm%GetMean(),& ! in - btran_eff, & ! in - vcmax_z, & ! out - jmax_z, & ! out - kp_z, & ! out - gs0, & ! out - gs1, & ! out - gs2 ) ! out - - - if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then - hydr_k_lwp = EDPftvarcon_inst%hydr_k_lwp(ft) - else - hydr_k_lwp = 1._r8 - end if - - call LeafLayerPhotosynthesis( & ! - par_abs, & ! in - ft, & ! in - vcmax_z, & ! in - jmax_z, & ! in - kp_z, & ! in - gs0, & ! in - gs1, & ! in - gs2, & ! in - bc_in(s)%t_veg_pa(ifp), & ! in - bc_in(s)%forc_pbot, & ! in - bc_in(s)%cair_pa(ifp), & ! in - bc_in(s)%oair_pa(ifp), & ! in - bc_in(s)%esat_tv_pa(ifp), & ! in - gb_mol, & ! in - bc_in(s)%eair_pa(ifp), & ! in - mm_kco2, & ! in - mm_ko2, & ! in - co2_cpoint, & ! in - lmr_z(iv,ft,cl), & ! in - ci_tol, & ! in - psn_ll, & ! out - gstoma_ll, & ! out - anet_ll, & ! out - c13disc_ll, & ! out - co2_inter_c_utest, & ! out (unit tests) - solve_iter) ! out performance tracking - - ! Average output quantities across sunlit and shaded leaves - ! Convert from molar to velocity (umol /m**2/s) to (m/s) - gstoma = gstoma + area_frac*(gstoma_ll / vmol_cf) - - - psn_z(iv,ft,cl) = psn_z(iv,ft,cl) + area_frac * psn_ll - anet_av_z(iv,ft,cl) = anet_av_z(iv,ft,cl) + area_frac * anet_ll - c13disc_z(iv,ft,cl) = c13disc_z(iv,ft,cl) + area_frac * c13disc_ll - - - end do do_sunsha - - ! Stomatal resistance of the leaf-layer - if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then - - rs_z(iv,ft,cl) = LeafHumidityStomaResis(leaf_psi, EDPftvarcon_inst%hydr_k_lwp(ft), & - bc_in(s)%t_veg_pa(ifp),bc_in(s)%cair_pa(ifp),bc_in(s)%forc_pbot, & - bc_in(s)%rb_pa(ifp), gstoma, ft, bc_in(s)%esat_tv_pa(ifp) ) - - else - rs_z(iv,ft,cl)= 1._r8/gstoma - end if - - rate_mask_z(iv,ft,cl) = .true. - - end if rate_mask_if - end do leaf_layer_loop - - - ! Zero cohort flux accumulators. - currentCohort%resp_m_tstep = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%ts_net_uptake = 0.0_r8 - currentCohort%c13disc_clm = 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 - - if(hlm_radiation_model.eq.norman_solver) then - - call ScaleLeafLayerFluxToCohort(nv, & !in - psn_z(1:nv,ft,cl), & !in - lmr_z(1:nv,ft,cl), & !in - rs_z(1:nv,ft,cl), & !in - currentPatch%elai_profile(cl,ft,1:nv), & !in - c13disc_z(1:nv,ft,cl), & !in - currentCohort%c_area, & !in - currentCohort%n, & !in - bc_in(s)%rb_pa(ifp), & !in - maintresp_reduction_factor, & !in - currentCohort%g_sb_laweight, & !out - currentCohort%gpp_tstep, & !out - currentCohort%rdark, & !out - currentCohort%c13disc_clm, & !out - cohort_eleaf_area) !out - - else - - - - call ScaleLeafLayerFluxToCohort(nv, & !in - psn_z(1:nv,ft,cl), & !in - lmr_z(1:nv,ft,cl), & !in - rs_z(1:nv,ft,cl), & !in - cohort_layer_elai(1:nv), & !in - c13disc_z(1:nv,ft,cl), & !in - currentCohort%c_area, & !in - currentCohort%n, & !in - bc_in(s)%rb_pa(ifp), & !in - maintresp_reduction_factor, & !in - currentCohort%g_sb_laweight, & !out - currentCohort%gpp_tstep, & !out - currentCohort%rdark, & !out - currentCohort%c13disc_clm, & !out - cohort_eleaf_area) !out - end if - - - ! Net Uptake does not need to be scaled, just transfer directly - currentCohort%ts_net_uptake(1:nv) = anet_av_z(1:nv,ft,cl) * umolC_to_kgC - - else - - ! In this case, the cohort had no leaves, - ! so no productivity,conductance, transpiration uptake - ! or dark respiration - cohort_eleaf_area = 0.0_r8 - currentCohort%gpp_tstep = 0.0_r8 - currentCohort%rdark = 0.0_r8 - currentCohort%g_sb_laweight = 0.0_r8 - currentCohort%ts_net_uptake(:) = 0.0_r8 + !write (fates_log(),*)'error, incorrect leaf respiration model specified' + !call endrun(msg=errMsg(sourcefile, __LINE__)) - end if canopy_mask_if + end select + ! Pre-process PAR absorbed per unit leaf area for different schemes + ! par_per_sunla = [W absorbed beam+diffuse radiation / m2 of sunlit leaves] + ! par_per_shala = [W absorbed diffuse radiation / m2 of shaded leaves] + ! fsun = [m2 of sunlit leaves / m2 of total leaves] + ! laisun: m2 of exposed leaf, per m2 of crown. If this is the lowest layer + ! for the pft/canopy group, than the m2 per crown is probably not + ! as large as the layer above. + ! ------------------------------------------------------------------ - ! ------------------------------------------------------------------ - ! Part VIII: Calculate maintenance respiration in the sapwood and - ! fine root pools. - ! ------------------------------------------------------------------ + if_radsolver: if(hlm_radiation_model.eq.norman_solver) then - ! Calculate the amount of nitrogen in the above and below ground - ! stem and root pools, used for maint resp - ! We are using the fine-root C:N ratio as an approximation for - ! the sapwood pools. - ! Units are in (kgN/plant) - ! ------------------------------------------------------------------ + laisun = patch%ed_laisun_z(cl,ft,iv) + laisha = patch%ed_laisha_z(cl,ft,iv) + par_per_sunla = patch%ed_parsun_z(cl,ft,iv) + par_per_shala = patch%ed_parsha_z(cl,ft,iv) + canopy_area = patch%canopy_area_profile(cl,ft,iv) + fsun = patch%f_sun(cl,ft,iv) - sapw_c = currentCohort%prt%GetState(sapw_organ, carbon12_element) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, carbon12_element) + else ! Two-stream - if (hlm_use_tree_damage .eq. itrue) then + if(cohort_layer_elai(iv) > nearzero .and. site%coszen>0._r8 ) then - ! Crown damage currenly only reduces the aboveground portion of - ! sapwood. Therefore we calculate the aboveground and the belowground portion - ! sapwood for use in stem respiration. - call GetCrownReduction(currentCohort%crowndamage, crown_reduction) + ! Since this is supposed to be cohort agnostic, lets + ! try to call the twostream function directly? + ! We should noot need to know cohort_elai here..? + call patch%twostr%GetAbsRad(cl,co_arr%twostr_col(ico),ipar,cohort_vaitop(iv), cohort_vaibot(iv), & + Rb_abs_el,Rd_abs_el,rd_abs_leaf,rb_abs_leaf,r_abs_stem_el,r_abs_snow_el,fsun,call_fail) + + !call FatesGetCohortAbsRad(patch, cohort, ipar, & + ! cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, & + ! rb_abs, rd_abs, rb_abs_leaf_el, rd_abs_leaf_el, leaf_sun_frac) + + !if(abs(rd_abs_leaf_el-rd_abs_leaf)>1.e-5) then + ! write(fates_log(),*)"Rds:",rd_abs_leaf_el,rd_abs_leaf + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + !end if + + !if(abs(rb_abs_leaf_el-rb_abs_leaf)>1.e-5) then + ! write(fates_log(),*)"Rbs:",rb_abs_leaf_el,rb_abs_leaf + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + !end if + + ! rd_abs_leaf: Watts of diffuse light absorbed by leaves over this + ! depth interval and ground footprint (m2) + ! rd_abs_leaf*fsun Watts of diffuse light absorbed by sunlit leaves + ! over this depth interval and ground footprint (m2) + ! rb_abs_leaf Watts of beam absorbed by sunlit leaves over this + ! depth interval and ground footprint (m2) + ! cohort_layer_elai*fsun Leaf area in sunlight within this interval and ground footprint + ! cohort_layer_elai*(1-fsun) Leaf area in shade within this interval and ground footprint + + laisun = (fsun*cohort_layer_elai(iv)) + laisha = ((1._r8 - fsun)*cohort_layer_elai(iv)) + if(fsun>nearzero) then + par_per_sunla = (rd_abs_leaf*fsun + rb_abs_leaf)! / laisun else - crown_reduction = 0.0_r8 + par_per_sunla = 0._r8 end if + par_per_shala = rd_abs_leaf*(1._r8-fsun) !/ laisha + canopy_area = 1._r8 !patch%canopy_area_profile(cl,ft,iv) + + else + + par_per_sunla = 0._r8 + par_per_shala = 0._r8 + laisun = 0.5_r8*cohort_layer_elai(iv) + laisha = 0.5_r8*cohort_layer_elai(iv) + canopy_area = 1._r8 !patch%canopy_area_profile(cl,ft,iv) + fsun = 0.5_r8 !avoid div0, should have no impact + + end if + + end if if_radsolver + + ! Perform photosynthesis calculations on sunlit and shaded leaves + ! --------------------------------------------------------------- + + ! Calculate leaf boundary layer conductance in molar form [umol/m2/s] + gb_mol = (1._r8/bc_in%rb_pa(ifp)) * vmol_cf + + ! 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(ft, & ! in + co_arr%vcmax25top(ico), & ! in + co_arr%jmax25top(ico), & ! in + co_arr%kp25top(ico), & ! in + nscaler, & ! in + bc_in%t_veg_pa(ifp), & ! in + bc_in%dayl_factor_pa(ifp), & ! in + patch%bprates, & ! in + btran_eff, & ! in + vcmax_z, & ! out + jmax_z, & ! out + kp_z, & ! out + gs0, & ! out + gs1, & ! out + gs2 ) ! out + + gstoma = 0._r8 + do_sunsha: do isunsha = 1,2 + + ! Determine absorbed PAR per square meter of leaf + ! If there is no leaf area perform a trivial solution + + if(isunsha == idirect) then + leaf_area = laisun*canopy_area + par_abs = ConvertPar(leaf_area, par_per_sunla) + area_frac = fsun + else + leaf_area = laisha*canopy_area + par_abs = ConvertPar(leaf_area, par_per_shala) + area_frac = 1._r8 - fsun + end if + + if( leaf_area < nearzero ) then + + ! Note: With no leaf area do not increment + ! any fluxes. Assume a nominal conductance + ! of maximum resistance + gstoma = gstoma + area_frac/rsmax0 + + cycle do_sunsha + end if + + if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then + hydr_k_lwp = EDPftvarcon_inst%hydr_k_lwp(ft) + else + hydr_k_lwp = 1._r8 + end if + + call LeafLayerPhotosynthesis( & ! + par_abs, & ! in + ft, & ! in + vcmax_z, & ! in + jmax_z, & ! in + kp_z, & ! in + gs0, & ! in + gs1, & ! in + gs2, & ! in + bc_in%t_veg_pa(ifp), & ! in + bc_in%forc_pbot, & ! in + bc_in%cair_pa(ifp), & ! in + bc_in%oair_pa(ifp), & ! in + bc_in%esat_tv_pa(ifp), & ! in + gb_mol, & ! in + bc_in%eair_pa(ifp), & ! in + mm_kco2, & ! in + mm_ko2, & ! in + co2_cpoint, & ! in + lmr_z(iv,uft,cl), & ! in + ci_tol, & ! in + psn_ll, & ! out + gstoma_ll, & ! out + anet_ll, & ! out + c13disc_ll, & ! out + co2_inter_c, & ! out (unit tests) + solve_iter) ! out performance tracking + + + ! Average output quantities across sunlit and shaded leaves + ! Convert from molar to velocity (umol /m**2/s) to (m/s) + gstoma = gstoma + area_frac*(gstoma_ll / vmol_cf) + + + psn_z(iv,uft,cl) = psn_z(iv,uft,cl) + area_frac * psn_ll + anet_av_z(iv,uft,cl) = anet_av_z(iv,uft,cl) + area_frac * anet_ll + c13disc_z(iv,uft,cl) = c13disc_z(iv,uft,cl) + area_frac * c13disc_ll + + + end do do_sunsha + + ! Stomatal resistance of the leaf-layer + if ( (hlm_use_planthydro.eq.itrue .and. EDPftvarcon_inst%hydr_k_lwp(ft)>nearzero) ) then + + rs_z(iv,uft,cl) = LeafHumidityStomaResis(leaf_psi, EDPftvarcon_inst%hydr_k_lwp(ft), & + bc_in%t_veg_pa(ifp),bc_in%cair_pa(ifp),bc_in%forc_pbot, & + bc_in%rb_pa(ifp), gstoma, ft, bc_in%esat_tv_pa(ifp) ) + + else + rs_z(iv,uft,cl)= 1._r8/gstoma + end if + + rate_mask_z(iv,uft,cl) = .true. + + end if rate_mask_if + end do leaf_layer_loop + + + ! Zero cohort flux accumulators. + co_arr%gpp_tstep(ico) = 0.0_r8 + co_arr%rdark_tstep(ico) = 0.0_r8 + co_arr%ts_net_uptake(:,ico) = 0.0_r8 + co_arr%c13disc_clm(ico) = 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 = co_arr%nv(ico) + + if(hlm_radiation_model.eq.norman_solver) then + + call ScaleLeafLayerFluxToCohort(nv, & !in + psn_z(1:nv,uft,cl), & !in + lmr_z(1:nv,uft,cl), & !in + rs_z(1:nv,uft,cl), & !in + patch%elai_profile(cl,ft,1:nv), & !in + c13disc_z(1:nv,uft,cl), & !in + co_arr%c_area(ico), & !in + co_arr%nplant(ico), & !in + bc_in%rb_pa(ifp), & !in + co_arr%mr_reduction_factor(ico), & !in + co_arr%g_sb_laweight(ico), & !out + co_arr%gpp_tstep(ico), & !out + co_arr%rdark_tstep(ico), & !out + co_arr%c13disc_clm(ico), & !out + cohort_eleaf_area) !out + + else + + call ScaleLeafLayerFluxToCohort(nv, & !in + psn_z(1:nv,uft,cl), & !in + lmr_z(1:nv,uft,cl), & !in + rs_z(1:nv,uft,cl), & !in + cohort_layer_elai(1:nv), & !in + c13disc_z(1:nv,uft,cl), & !in + co_arr%c_area(ico), & !in + co_arr%nplant(ico), & !in + bc_in%rb_pa(ifp), & !in + co_arr%mr_reduction_factor(ico), & !in + co_arr%g_sb_laweight(ico), & !out + co_arr%gpp_tstep(ico), & !out + co_arr%rdark_tstep(ico), & !out + co_arr%c13disc_clm(ico), & !out + cohort_eleaf_area) !out + end if + + + ! Net Uptake does not need to be scaled, just transfer directly + co_arr%ts_net_uptake(1:nv,ico) = anet_av_z(1:nv,uft,cl) * umolC_to_kgC + + ! Accumulate the combined conductance (stomatal+leaf boundary layer) + ! Note that cohort%g_sb_laweight is weighted by the leaf area + ! of each cohort and has units of [m/s] * [m2 leaf] + g_sb_leaves = g_sb_leaves + co_arr%g_sb_laweight(ico) - ! If crown reduction is zero, undamaged sapwood target will equal sapwood carbon - agb_frac = prt_params%allom_agb_frac(currentCohort%pft) - branch_frac = param_derived%branch_frac(currentCohort%pft) - sapw_c_undamaged = sapw_c / (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) - - ! Undamaged below ground portion - sapw_c_bgw = sapw_c_undamaged * (1.0_r8 - agb_frac) - - ! Damaged aboveground portion - sapw_c_agw = sapw_c - sapw_c_bgw - - - select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) - - live_stem_n = sapw_c_agw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - - live_croot_n = sapw_c_bgw * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - - fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) - - case(prt_cnp_flex_allom_hyp) - - live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & - currentCohort%prt%GetState(sapw_organ, nitrogen_element) - - live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & - currentCohort%prt%GetState(sapw_organ, nitrogen_element) - - - fnrt_n = currentCohort%prt%GetState(fnrt_organ, nitrogen_element) - - if (hlm_use_tree_damage .eq. itrue) then - - sapw_n = currentCohort%prt%GetState(sapw_organ, nitrogen_element) - - sapw_n_undamaged = sapw_n / & - (1.0_r8 - (agb_frac * branch_frac * crown_reduction)) - - sapw_n_bgw = sapw_n_undamaged * (1.0_r8 - agb_frac) - sapw_n_agw = sapw_n - sapw_n_bgw - - live_croot_n = sapw_n_bgw - - live_stem_n = sapw_n_agw - - end if - - ! If one wants to break coupling with dynamic N conentrations, - ! use the stoichiometry parameter - ! - ! live_stem_n = prt_params%allom_agb_frac(currentCohort%pft) * & - ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - ! live_croot_n = (1.0_r8-prt_params%allom_agb_frac(currentCohort%pft)) * & - ! sapw_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(sapw_organ)) - ! fnrt_n = fnrt_c * prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(fnrt_organ)) - - - case default - + ! Accumulate the total effective leaf area from all cohorts + ! in this patch. Normalize by canopy area outside the loop + patch_la = patch_la + cohort_eleaf_area - end select - !------------------------------------------------------------------------------ - ! 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 dealt with. - !------------------------------------------------------------------------------ + else - ! Live stem MR (kgC/plant/s) (above ground sapwood) - ! ------------------------------------------------------------------ - if ( int(woody(ft)) == itrue) then - tcwood = q10_mr**((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 * maintresp_nonleaf_baserate * tcwood * maintresp_reduction_factor - else - currentCohort%livestem_mr = 0._r8 - end if + ! In this case, the cohort had no leaves, + ! so no productivity,conductance, transpiration uptake + ! or dark respiration + cohort_eleaf_area = 0.0_r8 + co_arr%gpp_tstep(ico) = 0.0_r8 + co_arr%rdark_tstep(ico) = 0.0_r8 + co_arr%g_sb_laweight(ico) = 0.0_r8 + co_arr%ts_net_uptake(:,ico) = 0.0_r8 + end if canopy_mask_if + + enddo - ! Fine Root MR (kgC/plant/s) - ! and calculate the N fixation rate as a function of the fixation-specific root respiration - ! for now use dev_arbitrary_pft as scaling term between 0 and 1 as additional increment of root respiration used for N fixation - ! ------------------------------------------------------------------ - currentCohort%froot_mr = 0._r8 - currentCohort%sym_nfix_tstep = 0._r8 - ! n_fixation is integrated over the course of the day - ! this variable is zeroed at the end of the FATES dynamics sequence + ! Normalize canopy total conductance by the effective LAI + ! The value here was integrated over each cohort x leaf layer + ! and was weighted by m2 of effective leaf area for each layer - do j = 1,bc_in(s)%nlevsoil - tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) + if_any_lai: if(patch_la>nearzero) then + + ! Normalize the leaf-area weighted canopy conductance + ! The denominator is the total effective leaf area in the canopy, + ! units of [m/s]*[m2] / [m2] = [m/s] + + g_sb_leaves = g_sb_leaves / patch_la - fnrt_mr_layer = fnrt_n * maintresp_nonleaf_baserate * tcsoi * rootfr_ft(ft,j) * maintresp_reduction_factor + if_above_mincond: if( g_sb_leaves > (1._r8/rsmax0) ) then + + ! Combined mean leaf resistance is + ! the inverse of mean leaf conductance + r_sb_leaves = 1.0_r8/g_sb_leaves + + if (r_sb_leaves s/m (ideal gas conversion) [umol/m3] - currentCohort%sym_nfix_tstep = currentCohort%sym_nfix_tstep + nfix_layer + patch%c_stomata = vmol_cf / r_stomata + else !if_any_lai - enddo + ! But this will prevent it from using an unintialized value + bc_out%rssun_pa(ifp) = rsmax0 + bc_out%rssha_pa(ifp) = rsmax0 - ! Coarse Root MR (kgC/plant/s) (below ground sapwood) - ! ------------------------------------------------------------------ - if ( int(woody(ft)) == itrue) then - currentCohort%livecroot_mr = 0._r8 - do j = 1,bc_in(s)%nlevsoil - ! Soil temperature used to adjust base rate of MR - tcsoi = q10_mr**((bc_in(s)%t_soisno_sl(j)-tfrz - 20.0_r8)/10.0_r8) - currentCohort%livecroot_mr = currentCohort%livecroot_mr + & - live_croot_n * maintresp_nonleaf_baserate * tcsoi * & - rootfr_ft(ft,j) * maintresp_reduction_factor - enddo - else - currentCohort%livecroot_mr = 0._r8 - end if + ! This value is used for diagnostics, the molar form of conductance + ! is what is used in the field usually, so we track that form + patch%c_stomata = vmol_cf / rsmax0 + end if if_any_lai - ! ------------------------------------------------------------------ - ! Part IX: Perform some unit conversions (rate to integrated) and - ! calcualate some fluxes that are sums and nets of the base fluxes - ! ------------------------------------------------------------------ + ! This value is used for diagnostics, the molar form of conductance + ! is what is used in the field usually, so we track that form + patch%c_lblayer = vmol_cf / bc_in%rb_pa(ifp) - ! add on whole plant respiration values in kgC/indiv/s-1 - currentCohort%resp_m_tstep = currentCohort%livestem_mr + & - currentCohort%livecroot_mr + & - currentCohort%froot_mr + & - currentCohort%rdark - - ! no drought response right now.. something like: - ! resp_m_tstep = resp_m_tstep * (1.0_r8 - currentPatch%btran_ft(currentCohort%pft) * & - ! EDPftvarcon_inst%resp_drought_response(ft)) - - ! convert from kgC/indiv/s to kgC/indiv/timestep - currentCohort%resp_m_tstep = currentCohort%resp_m_tstep * dtime - currentCohort%gpp_tstep = currentCohort%gpp_tstep * dtime - currentCohort%ts_net_uptake = currentCohort%ts_net_uptake * dtime - - ! save as a diagnostic the un-throttled maintenance respiration to be able to know how strong this is - currentCohort%resp_m_unreduced = currentCohort%resp_m_tstep / maintresp_reduction_factor - - ! Accumulate the combined conductance (stomatal+leaf boundary layer) - ! Note that currentCohort%g_sb_laweight is weighted by the leaf area - ! of each cohort and has units of [m/s] * [m2 leaf] - - g_sb_leaves = g_sb_leaves + currentCohort%g_sb_laweight - - ! Accumulate the total effective leaf area from all cohorts - ! in this patch. Normalize by canopy area outside the loop - patch_la = patch_la + cohort_eleaf_area - - currentCohort => currentCohort%shorter - enddo do_cohort_drive + end block + end associate + + return + end subroutine FatesCondPhotoPatch - end if if_any_cohorts + ! ============================================================================= - ! Normalize canopy total conductance by the effective LAI - ! The value here was integrated over each cohort x leaf layer - ! and was weighted by m2 of effective leaf area for each layer - - if_any_lai: if(patch_la>nearzero) then + subroutine FatesPlantRespPatch(ifp,site,bc_in,dtime) - ! Normalize the leaf-area weighted canopy conductance - ! The denominator is the total effective leaf area in the canopy, - ! units of [m/s]*[m2] / [m2] = [m/s] - - g_sb_leaves = g_sb_leaves / patch_la - - if_above_mincond: if( g_sb_leaves > (1._r8/rsmax0) ) then - - ! Combined mean leaf resistance is - ! the inverse of mean leaf conductance - r_sb_leaves = 1.0_r8/g_sb_leaves - - if (r_sb_leaves s/m (ideal gas conversion) [umol/m3] + ! ----------------------------------------------------------------------------------- + ! !DESCRIPTION: + ! Plant respiration (non-dark) from various organs. + ! This is separated from photosynthesis and conductance because those + ! calculations may need to be called multiple times, as they are + ! often called inside a land-energy balance solver + ! ----------------------------------------------------------------------------------- - currentPatch%c_stomata = vmol_cf / r_stomata - - else !if_any_lai - - ! But this will prevent it from using an unintialized value - bc_out(s)%rssun_pa(ifp) = rsmax0 - bc_out(s)%rssha_pa(ifp) = rsmax0 - - ! This value is used for diagnostics, the molar form of conductance - ! is what is used in the field usually, so we track that form - currentPatch%c_stomata = vmol_cf / rsmax0 - - end if if_any_lai + integer :: ifp + + type(ed_site_type),intent(inout) :: site + type(bc_in_type),intent(in) :: bc_in + real(r8),intent(in) :: dtime + + ! Locals + type (fates_patch_type), pointer :: patch + real(r8) :: fnrt_mr_nfix_layer ! fineroot maintenance respiration + ! specifically for symbiotic fixation [kgC/plant/layer/s] + real(r8) :: nfix_layer ! Nitrogen fixed in each layer this timestep [kgN/plant/layer/timestep] + + real(r8) :: tcsoi(bc_in%nlevsoil) ! Temperature response function for root respiration. + real(r8) :: tcwood ! Temperature response function for wood + real(r8) :: fnrt_mr_layer ! fine root maintenance respiation per layer [kgC/plant/s] + integer :: ft, j ! indices: pft, soil layer + real(r8) :: fnrtfrac_ftz(numpft,bc_in%nlevsoil) + integer :: ico + + patch => site%pa_vec(ifp)%p - ! This value is used for diagnostics, the molar form of conductance - ! is what is used in the field usually, so we track that form - currentPatch%c_lblayer = vmol_cf / bc_in(s)%rb_pa(ifp) - - end if if_filter2 - - end if if_notbare + associate(coarr => patch%coarrays) + + ! ----------------------------------------------------------------------------------- + ! Calculate the amount of nitrogen in the above and below ground stem and root + ! pools, used for maint resp. We are using the fine-root C:N ratio as an + ! approximation for the sapwood pools. Units are in (kgN/plant) + ! ----------------------------------------------------------------------------------- - currentPatch => currentPatch%younger - end do + do ft = 1,numpft + call set_root_fraction(fnrtfrac_ftz(ft,:), ft, & + bc_in%zi_sisl, & + bc_in%max_rooting_depth_index_col) + end do - deallocate(rootfr_ft) + ! Temperature correction for wood on this patch + tcwood = FastQ10(bc_in%t_veg_pa(patch%patchno),20._r8,q10_mr,log_q10_mr_div10) + do j = 1,bc_in%nlevsoil + tcsoi(j) = FastQ10(bc_in%t_soisno_sl(j),20._r8,q10_mr,log_q10_mr_div10) + end do - end do !site loop - end associate - end subroutine FatesPlantRespPhotosynthDrive + do_cohort: do ico = 1,co_arr%ncohorts + + ft = coarr%pft(ico) + + ! Live stem MR (kgC/plant/s) (above ground sapwood) + ! ------------------------------------------------------------------ + ! kgC/s = kgN * kgC/kgN/s + + coarr%livestem_mr(ico) = coarr%live_stem_n(ico) * maintresp_nonleaf_baserate * & + tcwood * patch%coarrays%mr_reduction_factor(ico) + + ! Fine Root MR (kgC/plant/s) + ! and calculate the N fixation rate as a function of the fixation-specific + ! root respiration for now use dev_arbitrary_pft as scaling term between 0 + ! and 1 as additional increment of root respiration used for N fixation + ! ------------------------------------------------------------------ + coarr%froot_mr(ico) = 0._r8 + coarr%sym_nfix_tstep(ico) = 0._r8 + + ! n_fixation is integrated over the course of the day + ! this variable is zeroed at the end of the FATES dynamics sequence + + do j = 1,bc_in%nlevsoil + + fnrt_mr_layer = coarr%fnrt_n(ico) * maintresp_nonleaf_baserate * tcsoi(j) * & + fnrtfrac_ftz(ft,j) * coarr%mr_reduction_factor(ico) + + ! calculate the cost of carbon for N fixation in each soil layer + ! and calculate N fixation rate based on that [kgC / kgN] + + call RootLayerNFixation(bc_in%t_soisno_sl(j),ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,nfix_layer) + + coarr%froot_mr(ico) = coarr%froot_mr(ico) + fnrt_mr_nfix_layer + fnrt_mr_layer + + coarr%sym_nfix_tstep(ico) = coarr%sym_nfix_tstep(ico) + nfix_layer + + enddo + + ! Coarse Root MR (kgC/plant/s) (below ground sapwood) + ! ------------------------------------------------------------------ + if ( int(prt_params%woody(ft)) == itrue) then + coarr%livecroot_mr(ico) = 0._r8 + do j = 1,bc_in%nlevsoil + ! Soil temperature used to adjust base rate of MR + coarr%livecroot_mr(ico) = coarr%livecroot_mr(ico) + & + coarr%live_croot_n(ico) * maintresp_nonleaf_baserate * tcsoi(j) * & + fnrtfrac_ftz(ft,j) * coarr%mr_reduction_factor(ico) + enddo + else + coarr%livecroot_mr(ico) = 0._r8 + end if + + ! ------------------------------------------------------------------ + ! Part IX: Perform some unit conversions (rate to integrated) and + ! calcualate some fluxes that are sums and nets of the base fluxes + ! ------------------------------------------------------------------ + + ! add on whole plant respiration values in kgC/indiv/s-1 + coarr%resp_m_tstep(ico) = coarr%livestem_mr(ico) + & + coarr%livecroot_mr(ico) + & + coarr%froot_mr(ico) + & + coarr%rdark_tstep(ico) + + ! no drought response right now.. something like: + ! resp_m_tstep = resp_m_tstep * (1.0_r8 - patch%btran_ft(cohort%pft) * & + ! EDPftvarcon_inst%resp_drought_response(ft)) + + ! convert from kgC/indiv/s to kgC/indiv/timestep + coarr%resp_m_tstep(ico) = coarr%resp_m_tstep(ico) * dtime + coarr%gpp_tstep(ico) = coarr%gpp_tstep(ico) * dtime + coarr%ts_net_uptake(:,ico) = coarr%ts_net_uptake(:,ico) * dtime + + ! save as a diagnostic the un-throttled maintenance respiration to be able to know how strong this is + coarr%resp_m_unreduced(ico) = coarr%resp_m_tstep(ico) / coarr%mr_reduction_factor(ico) + + enddo do_cohort_drive + + end subroutine FatesPlantRespPatch ! =========================================================================================== - subroutine RootLayerNFixation(t_soil,ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,nfix_layer) @@ -1201,7 +1086,7 @@ subroutine RootLayerNFixation(t_soil,ft,dtime,fnrt_mr_layer,fnrt_mr_nfix_layer,n ! Amount of nitrogen fixed in this layer [kgC/plant/layer/tstep]/[kgC/kgN] = [kgN/plant/layer/tstep] nfix_layer = c_spent_nfix / c_cost_nfix - + return end subroutine RootLayerNFixation @@ -1254,49 +1139,15 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in real(r8) :: sum_weight ! sum of weight for unpacking d13c flux (c13disc_z) from ! (canopy_layer, pft, leaf_layer) matrix to cohort (c13disc_clm) - ! 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 integer :: il ! leaf layer index real(r8) :: cohort_layer_eleaf_area ! the effective leaf area of the cohort's current layer [m2] - cohort_eleaf_area = 0.0_r8 - g_sb_laweight = 0.0_r8 - gpp = 0.0_r8 - rdark = 0.0_r8 - - do il = 1, nv ! Loop over the leaf layers this cohort participates in - - - ! Cohort's total effective leaf area in this layer [m2] - ! leaf area index of the layer [m2/m2 ground] * [m2 ground] - ! elai_llz is the LAI for the whole PFT. Multiplying this by the ground - ! area this cohort contributes, give the cohort's portion of the leaf - ! area in this layer - cohort_layer_eleaf_area = elai_llz(il) * c_area - - ! Increment the cohort's total effective leaf area [m2] - cohort_eleaf_area = cohort_eleaf_area + cohort_layer_eleaf_area - - ! Leaf conductance (stomatal and boundary layer) - ! This should be the weighted average over the leaf surfaces. - ! Since this is relevant to the stomata, its weighting should be based - ! on total leaf area, and not really footprint area - ! [m/s] * [m2 cohort's leaf layer] - g_sb_laweight = g_sb_laweight + 1.0_r8/(rs_llz(il)+rb) * cohort_layer_eleaf_area - - ! GPP [umolC/m2leaf/s] * [m2 leaf ] -> [umolC/s] - gpp = gpp + psn_llz(il) * cohort_layer_eleaf_area - - ! Dark respiration - ! [umolC/m2leaf/s] * [m2 leaf] - rdark = rdark + lmr_llz(il) * cohort_layer_eleaf_area - - end do - - - + gpp = sum(psn_llz(1:nv) * elai_llz(1:nv))*c_area + g_sb_laweight = sum(1.0_r8/(rs_llz(1:nv)+rb) * elai_llz(1:nv))*c_area + rdark = sum(lmr_llz(1:nv) * elai_llz(1:nv))*c_area + cohort_eleaf_area = sum(elai_llz(1:nv))*c_area + if (nv > 1) then ! cohort%c13disc_clm as weighted mean of d13c flux at all related leave layers sum_weight = sum(psn_llz(1:nv-1) * elai_llz(1:nv-1)) @@ -1324,90 +1175,11 @@ subroutine ScaleLeafLayerFluxToCohort(nv, & ! in rdark = rdark * umolC_to_kgC * maintresp_reduction_factor / nplant gpp = gpp * umolC_to_kgC / nplant - 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 873 ', nv - write(fates_log(),*) 'EDPhoto 874 ', cohort_eleaf_area - endif - return end subroutine ScaleLeafLayerFluxToCohort ! ===================================================================================== - subroutine UpdateCanopyNCanNRadPresent(currentPatch) - - ! --------------------------------------------------------------------------------- - ! This subroutine calculates two patch level quanities: - ! currentPatch%ncan and - ! currentPatch%canopy_mask - ! - ! 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%canopy_mask(:,:) has the same dimensions, is binary, and - ! indicates whether or not leaf layers are present (by evaluating the canopy area - ! profile). - ! --------------------------------------------------------------------------------- - - ! Arguments - type(fates_patch_type), target :: currentPatch - type(fates_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%nleaf(:,:) = 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%nleaf(currentCohort%canopy_layer,currentCohort%pft) = & - max(currentPatch%nleaf(currentCohort%canopy_layer,currentCohort%pft), & - currentCohort%NV) - - currentCohort => currentCohort%shorter - - enddo !cohort - - ! NRAD = NCAN ... - currentPatch%nrad = currentPatch%nleaf - - ! Now loop through and identify which layer and pft combo has scattering elements - do cl = 1,nclmax - do ft = 1,numpft - currentPatch%canopy_mask(cl,ft) = 0 - do iv = 1, currentPatch%nrad(cl,ft); - if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then - currentPatch%canopy_mask(cl,ft) = 1 - end if - end do !iv - enddo !ft - enddo !cl - - return - end subroutine UpdateCanopyNCanNRadPresent - - ! ===================================================================================== - - real(r8) function ConvertPar(leaf_area, par_wm2) result(par_umolm2s) ! @@ -1432,4 +1204,9 @@ real(r8) function ConvertPar(leaf_area, par_wm2) result(par_umolm2s) end function ConvertPar + ! ================================================================================= + + + + end module FATESPlantRespPhotosynthMod diff --git a/biogeophys/LeafBiophysicsMod.F90 b/biogeophys/LeafBiophysicsMod.F90 index c238eebf16..3144c4ceb2 100644 --- a/biogeophys/LeafBiophysicsMod.F90 +++ b/biogeophys/LeafBiophysicsMod.F90 @@ -64,11 +64,16 @@ module LeafBiophysicsMod public :: CiMinMax public :: CiFunc public :: CiBisection + public :: UpdateSlowBiophysicalRates + public :: UpdateFastBiophysicalRates + public :: FastQ10 character(len=*), parameter, private :: sourcefile = & __FILE__ + logical, parameter :: do_b4b = .false. + character(len=1024) :: warn_msg ! for defining a warning message !------------------------------------------------------------------------------------- @@ -116,14 +121,16 @@ module LeafBiophysicsMod integer, public, parameter :: JohnsonBerry2021 = 2 ! Constants defining the photosynthesis temperature acclimation model - integer, parameter :: photosynth_acclim_model_none = 0 - integer, parameter :: photosynth_acclim_model_kumarathunge_etal_2019 = 1 + integer, public, parameter :: photosynth_acclim_model_none = 0 + integer, public, parameter :: photosynth_acclim_model_kumarathunge_etal_2019 = 1 ! Rdark constants from Atkin et al., 2017 https://doi.org/10.1007/978-3-319-68703-2_6 ! and Heskel et al., 2016 https://doi.org/10.1073/pnas.1520282113 real(r8), parameter :: lmr_b = 0.1012_r8 ! (degrees C**-1) real(r8), parameter :: lmr_c = -0.0005_r8 ! (degrees C**-2) real(r8), parameter :: lmr_TrefC = 25._r8 ! (degrees C) + + real(r8), parameter :: t_ref_298K = tfrz + 25._r8 ! (25C in K) ! These two are public for error checking during parameter read-in real(r8), parameter, public :: lmr_r_1 = 0.2061_r8 ! (umol CO2/m**2/s / (gN/(m2 leaf))) @@ -136,7 +143,6 @@ module LeafBiophysicsMod ! empirical curvature parameter for electron transport rate real(r8),parameter :: theta_psii = 0.7_r8 - ! curvature parameter for quadratic smoothing of C4 gross assimilation real(r8),parameter :: theta_ip_c4 = 0.95_r8 !0.95 is from Collatz 91, 0.999 was api 36 real(r8),parameter :: theta_cj_c4 = 0.98_r8 !0.98 from Collatz 91, 0.099 was api 36 @@ -156,8 +162,41 @@ module LeafBiophysicsMod ! For plants with no leaves, a miniscule amount of conductance ! can happen through the stems, at a partial rate of cuticular conductance real(r8),parameter :: stem_cuticle_loss_frac = 0.1_r8 + + type, public :: bprate_type + + ! This holds various biophysical rates used in photosynthesis + ! These values change at slow timescales when using Kumarathunge et al. + ! The temperature response and inhibition functions require non-trivial + ! mathematical operations at high frequencies, so we pre-calculate + ! as much as we can before operating on leaf and organ temperatures. + ! For leaves, many of these pre-calculated rates are either constant + ! or slowly acclimating following Kumarathunge et al. + + real(r8),allocatable :: vr_a(:) ! term a of optimized temp response ft1_f + real(r8),allocatable :: vr_b(:) ! term b of optimized temp response ft1_f + real(r8),allocatable :: jr_a(:) ! term a of optimized temp response ft1_f + real(r8),allocatable :: jr_b(:) ! term b of optimized temp response ft1_f + real(r8),allocatable :: vi_a(:) ! term a of optimized temp inhibition fth1_f + real(r8),allocatable :: vi_b(:) ! term b of optimized temp inhibition fth1_f + real(r8),allocatable :: vi_c(:) ! term c of optimized temp inhibition fth1_f + real(r8),allocatable :: ji_a(:) ! term a of optimized temp inhibition fth1_f + real(r8),allocatable :: ji_b(:) ! term b of optimized temp inhibition fth1_f + real(r8),allocatable :: ji_c(:) ! term c of optimized temp inhibition fth1_f + real(r8),allocatable :: jvr(:) ! ratio of Jmax25 / Vcmax25 ! Kumarathunge et al. + real(r8),allocatable :: vcmax_tscaler(:) ! Temperature scaling functions on vcmax + real(r8),allocatable :: jmax_tscaler(:) ! Temperature scaling functions on jmax + real(r8),allocatable :: lmr_tscaler(:) ! Temperature response scaler on LMR + + contains + procedure :: BPRateInitAlloc + procedure :: BPRateDealloc + end type bprate_type + + + ! The stomatal slope can either be scaled by btran or not. FATES had ! a precedent of using this into 2024, but discussions here: https://github.com/NGEET/fates/issues/719 ! suggest we should try other hypotheses @@ -223,6 +262,14 @@ module LeafBiophysicsMod ! 1: btran scales only vcmax ! 2: btran scales both vcmax and jmax real(r8),allocatable :: fnps(:) ! fraction of light absorbed by non-photosynthetic pigments + + + + ! Pre-calculations for improved performance + + real(r8),allocatable :: vcmaxc_25(:) + real(r8),allocatable :: jmaxc_25(:) + ! ------------------------------------------------------------------------------------- ! Note the omission of several parameter constants: ! @@ -240,8 +287,6 @@ module LeafBiophysicsMod type(leafbiophys_params_type),public :: lb_params - - ! A possible sequence of calls for leaf biophysics is as follows: ! 1) determine any gas quantities or parameters that are derived and ! are applicable to a super-set of leaf-layers (like MM, and compensation points) @@ -254,6 +299,57 @@ module LeafBiophysicsMod contains + subroutine BPRateInitAlloc(this,numpft) + + ! Class helper routine that allocates the type holding + ! biophyical rates + + class(bprate_type) :: this + integer,intent(in) :: numpft + + allocate(this%vr_a(numpft));this%vr_a(:)=fates_unset_r8 + allocate(this%vr_b(numpft));this%vr_b(:)=fates_unset_r8 + allocate(this%jr_a(numpft));this%jr_a(:)=fates_unset_r8 + allocate(this%jr_b(numpft));this%jr_b(:)=fates_unset_r8 + allocate(this%vi_a(numpft));this%vi_a(:)=fates_unset_r8 + allocate(this%vi_b(numpft));this%vi_b(:)=fates_unset_r8 + allocate(this%vi_c(numpft));this%vi_c(:)=fates_unset_r8 + allocate(this%ji_a(numpft));this%ji_a(:)=fates_unset_r8 + allocate(this%ji_b(numpft));this%ji_b(:)=fates_unset_r8 + allocate(this%ji_c(numpft));this%ji_c(:)=fates_unset_r8 + allocate(this%jvr(numpft));this%jvr(:)=fates_unset_r8 + allocate(this%vcmax_tscaler(numpft));this%vcmax_tscaler(:)=fates_unset_r8 + allocate(this%jmax_tscaler(numpft));this%jmax_tscaler(:)=fates_unset_r8 + allocate(this%lmr_tscaler(numpft));this%lmr_tscaler(:)=fates_unset_r8 + return + end subroutine BPRateInitAlloc + + ! ===================================================================================== + + subroutine BPRateDealloc(this) + + ! Class helper routine that deallocates biophysical rates + + class(bprate_type) :: this + deallocate(this%vr_a, & + this%vr_b, & + this%jr_a, & + this%jr_b, & + this%vi_a, & + this%vi_b, & + this%vi_c, & + this%ji_a, & + this%ji_b, & + this%ji_c, & + this%jvr, & + this%vcmax_tscaler, & + this%jmax_tscaler, & + this%lmr_tscaler) + return + end subroutine BPRateDealloc + + ! ===================================================================================== + subroutine StomatalCondMedlyn(anet,veg_esat,can_vpress,gs0,gs1,gs2,leaf_co2_ppress,can_press,gb,gs) ! ----------------------------------------------------------------------------------- @@ -320,13 +416,15 @@ subroutine StomatalCondMedlyn(anet,veg_esat,can_vpress,gs0,gs1,gs2,leaf_co2_ppre vpd = max((veg_esat - can_vpress), min_vpd_pa) * kpa_per_pa term = gs2 * h2o_co2_stoma_diffuse_ratio * anet / (leaf_co2_ppress / can_press) - aquad = 1.0_r8 + aquad = 1._r8 bquad = -(2.0 * (gs0 + term) + (gs1 * term)**2 / (gb * vpd )) cquad = gs0*gs0 + (2.0*gs0 + term * & (1.0 - gs1*gs1 / vpd)) * term - - call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) - gs = max(r1,r2) + + gs = QuadraticSolveSmooth_a1_max(bquad, cquad) + + !call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + !gs = max(r1,r2) if(debug)then if(err)then @@ -362,22 +460,22 @@ subroutine StomatalCondBallBerry(a_gs,a_net,veg_esat,can_vpress,gs0,gs1,leaf_co2 ! Input - real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) - real(r8), intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) - real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol /m**2/s) - real(r8), intent(in) :: can_vpress ! vapor pressure of the canopy air (Pa) - real(r8), intent(in) :: gs0,gs1 ! Stomatal intercept and slope (umol H2O/m**2/s) - real(r8), intent(in) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) - real(r8), intent(in) :: a_gs ! The assimilation (a) for calculating conductance (gs) - ! is either = to anet or agross - real(r8), intent(in) :: a_net ! Net assimilation rate of co2 (umol/m2/s) - ! Output - real(r8), intent(out) :: gs ! leaf stomatal conductance (umol H2O/m**2/s) - - ! locals - real(r8) :: ceair ! constrained vapor pressure (Pa) - real(r8) :: aquad,bquad,cquad ! quadradic solve terms - real(r8) :: r1,r2 ! quadradic solve roots + real(r8), intent(in) :: veg_esat ! saturation vapor pressure at veg_tempk (Pa) + real(r8), intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) + real(r8), intent(in) :: gb ! leaf boundary layer conductance (umol /m**2/s) + real(r8), intent(in) :: can_vpress ! vapor pressure of the canopy air (Pa) + real(r8), intent(in) :: gs0,gs1 ! Stomatal intercept and slope (umol H2O/m**2/s) + real(r8), intent(in) :: leaf_co2_ppress ! CO2 partial pressure at leaf surface (Pa) + real(r8), intent(in) :: a_gs ! The assimilation (a) for calculating conductance (gs) + ! is either = to anet or agross + real(r8), intent(in) :: a_net ! Net assimilation rate of co2 (umol/m2/s) + ! Output + real(r8), intent(out) :: gs ! leaf stomatal conductance (umol H2O/m**2/s) + + ! locals + real(r8) :: ceair ! constrained vapor pressure (Pa) + real(r8) :: aquad,bquad,cquad ! quadradic solve terms + real(r8) :: r1,r2 ! quadradic solve roots logical :: err @@ -403,16 +501,8 @@ subroutine StomatalCondBallBerry(a_gs,a_net,veg_esat,can_vpress,gs0,gs1,leaf_co2 else cquad = -gb*(leaf_co2_ppress * gs0 + gs1 * a_gs* can_press * ceair/ veg_esat ) end if - - call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) - if(debug)then - if(err)then - write(fates_log(),*) "bbquadfail:",a_net,a_gs,veg_esat,can_vpress,gs0,gs1,leaf_co2_ppress,can_press - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if - + call QuadraticSolveSmooth(aquad,bquad,cquad,r1,r2) gs = max(r1,r2) return @@ -420,15 +510,15 @@ end subroutine StomatalCondBallBerry ! ===================================================================================== - function AgrossRubiscoC3(vcmax,ci,can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) result(ac) + elemental function AgrossRubiscoC3(vcmax,ci,can_o2_ppress,co2_cpoint,mm_kco2,mm_ko2) result(ac) ! Input - real(r8) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) - real(r8) :: ci ! intracellular leaf CO2 (Pa) - real(r8) :: co2_cpoint ! CO2 compensation point (Pa) - real(r8) :: can_o2_ppress ! Partial pressure of O2 near the leaf surface (Pa) - real(r8) :: mm_kco2 ! Michaelis-Menten constant for CO2 (Pa) - real(r8) :: mm_ko2 ! Michaelis-Menten constant for O2 (Pa) + real(r8),intent(in) :: vcmax ! maximum rate of carboxylation (umol co2/m**2/s) + real(r8),intent(in) :: ci ! intracellular leaf CO2 (Pa) + real(r8),intent(in) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8),intent(in) :: can_o2_ppress ! Partial pressure of O2 near the leaf surface (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) ! Output real(r8) :: ac ! Rubisco-limited gross photosynthesis (umol CO2/m**2/s) @@ -440,13 +530,13 @@ end function AgrossRubiscoC3 ! ===================================================================================== - function GetJe(par_abs,jmax,fnps) result (je) + elemental function GetJe(par_abs,jmax,fnps) result (je) ! Input - real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] - real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments - real(r8) :: je ! electron transport rate (umol electrons/m**2/s) + real(r8),intent(in) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] + real(r8),intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8),intent(in) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments + real(r8) :: je ! electron transport rate (umol electrons/m**2/s) select case(lb_params%electron_transport_model) case (FvCB1980) @@ -456,23 +546,21 @@ function GetJe(par_abs,jmax,fnps) result (je) case (JohnsonBerry2021) je = GetJe_JB(par_abs,jmax,fnps) - case default - write (fates_log(),*)'error, incorrect leaf electron transport model specified:',lb_params%electron_transport_model - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select + end select end function GetJe !====================================================================================== - function GetJe_FvCB(par_abs,jmax,fnps) result(je) + elemental function GetJe_FvCB(par_abs,jmax,fnps) result(je) ! Input - real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] - real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments + real(r8),intent(in) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] + real(r8),intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8),intent(in) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments real(r8) :: je ! electron transport rate (umol electrons/m**2/s) + real(r8) :: aquad,bquad,cquad ! terms for quadratic equations real(r8) :: r1,r2 ! roots of quadratic equation real(r8) :: jpar ! absorbed photons in photocenters as an electron @@ -490,38 +578,41 @@ function GetJe_FvCB(par_abs,jmax,fnps) result(je) aquad = theta_psii bquad = -(jpar + jmax) cquad = jpar * jmax - call QuadraticRoots(aquad, bquad, cquad, r1, r2, err) - if(debug)then - if(err)then - write(fates_log(),*) "jequadfail:",par_abs,jpar,jmax - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + if(do_b4b) then + call QuadraticRoots(aquad, bquad, cquad, r1, r2, err) + !if(debug)then + ! if(err)then + ! write(fates_log(),*) "jequadfail:",par_abs,jpar,jmax + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + ! end if + !end if + je = min(r1,r2) + else + je = QuadraticSolveSmooth_theta_psii(bquad,cquad) end if - - je = min(r1,r2) - + end function GetJe_FvCB ! ===================================================================================== - function GetJe_JB(par_abs,jmax,fnps) result(je) + elemental function GetJe_JB(par_abs,jmax,fnps) result(je) ! Input - real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] - real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments - real(r8) :: je ! electron transport rate (umol electrons/m**2/s) - real(r8) :: phi ! maximum quantum yield (mol electrons/mol photons) - real(r8) :: cb6fmax ! maximum activity of the cytochrome b6f complex - ! (umol electrons/m**2/s) - ! referred to as vqmax in Lamour et al. - real(r8) :: Qsat ! Saturating irradiance - assumed to be a constant (umol/m**2/s) - ! Here we asssume abosorbed irradiance - real(r8) :: jsat ! Electron transport rate estimated by the FvCB model for a - ! given Jmax at Qsat - - Qsat = 1530.0_r8 + real(r8),intent(in) :: par_abs ! Absorbed PAR per leaf area [umol photons/m**2/s] + real(r8),intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8),intent(in) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments + real(r8) :: je ! electron transport rate (umol electrons/m**2/s) + + real(r8) :: phi ! maximum quantum yield (mol electrons/mol photons) + real(r8) :: cb6fmax ! maximum activity of the cytochrome b6f complex + ! (umol electrons/m**2/s) + ! referred to as vqmax in Lamour et al. + ! Here we asssume abosorbed irradiance + real(r8) :: jsat ! Electron transport rate estimated by the FvCB model for a + ! given Jmax at Qsat + real(r8),parameter :: Qsat = 1530.0_r8 ! Saturating irradiance - assumed to be + ! a constant (umol/m**2/s) phi = (1.0_r8 - fnps) * photon_to_e @@ -540,14 +631,14 @@ end function GetJe_JB ! ===================================================================================== - function AgrossRuBPC3(par_abs,jmax,fnps,ci,co2_cpoint) result(aj) + elemental function AgrossRuBPC3(par_abs,jmax,fnps,ci,co2_cpoint) result(aj) ! Input - real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] - real(r8) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) - real(r8) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments - real(r8) :: ci ! intracellular leaf CO2 (Pa) - real(r8) :: co2_cpoint ! CO2 compensation point (Pa) + real(r8),intent(in) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] + real(r8),intent(in) :: jmax ! maximum electron transport rate (umol electrons/m**2/s) + real(r8),intent(in) :: fnps ! Fraction of light absorbed by non-photosynthetic pigments + real(r8),intent(in) :: ci ! intracellular leaf CO2 (Pa) + real(r8),intent(in) :: co2_cpoint ! CO2 compensation point (Pa) ! Output real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) @@ -566,13 +657,13 @@ end function AgrossRuBPC3 ! ======================================================================================= - function AgrossRuBPC4(par_abs) result(aj) + elemental function AgrossRuBPC4(par_abs) result(aj) - real(r8) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] + real(r8),intent(in) :: par_abs ! Absorbed PAR per leaf area [umol photons/m2leaf/s ] real(r8) :: aj ! RuBP-limited gross photosynthesis (umol CO2/m**2/s) ! quantum efficiency, used only for C4 (mol CO2 / mol photons) - real(r8),parameter :: c4_quant_eff = 0.05_r8 + real(r8), parameter :: c4_quant_eff = 0.05_r8 aj = c4_quant_eff*par_abs @@ -580,11 +671,11 @@ end function AgrossRuBPC4 ! ======================================================================================= - function AgrossPEPC4(ci,kp,can_press) result(ap) + elemental function AgrossPEPC4(ci,kp,can_press) result(ap) - real(r8) :: ci ! intracellular leaf CO2 (Pa) - real(r8) :: kp ! initial co2 response slope - real(r8) :: can_press ! Air pressure near the surface of the leaf (Pa) + real(r8),intent(in) :: ci ! intracellular leaf CO2 (Pa) + real(r8),intent(in) :: kp ! initial co2 response slope + real(r8),intent(in) :: can_press ! Air pressure near the surface of the leaf (Pa) real(r8) :: ap ! PEP limited gross assimilation rate (umol co2/m2/s) ap = kp * max(ci, 0._r8) / can_press @@ -668,6 +759,8 @@ subroutine CiMinMax(ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, & ai = min(r1,r2) + + a = rmin*can_press aquad = theta_ip_c4/a**2.0_r8 + kp/(can_press*a) @@ -880,7 +973,7 @@ function CiFromAnetDiffGrad(a,b,c,d,e,f,g) result(ci) ! This function is called to find endpoints, where conductance is maximized ! and minimized, to perform a binary search - real(r8) :: a,b,c,d,e,f,g ! compound terms to solve the coupled Anet + real(r8),intent(in) :: a,b,c,d,e,f,g ! compound terms to solve the coupled Anet ! and diffusive flux gradient equations real(r8) :: ci ! intracellular co2 [Pa] real(r8) :: r1,r2 ! roots for quadratic @@ -990,30 +1083,30 @@ subroutine CiFunc(ci, & aquad = theta_cj_c4 bquad = -(ac + aj) cquad = ac * aj - call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) - - if(debug)then - if(err)then - write(fates_log(),*) "c41quadfail:",par_abs,ci,kp,can_press,vcmax - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if - - ai = min(r1,r2) - + !call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + !if(debug)then + ! if(err)then + ! write(fates_log(),*) "c41quadfail:",par_abs,ci,kp,can_press,vcmax + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + ! end if + !end if + !ai = min(r1,r2) + + ai = QuadraticSolveSmooth_theta_cj(bquad, cquad) + aquad = theta_ip_c4 bquad = -(ai + ap) cquad = ai * ap - call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) - - if(debug)then - if(err)then - write(fates_log(),*) "c42quadfail:",par_abs,ci,kp,can_press,vcmax - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if + !call QuadraticRoots(aquad, bquad, cquad, r1, r2,err) + !if(debug)then + ! if(err)then + ! write(fates_log(),*) "c42quadfail:",par_abs,ci,kp,can_press,vcmax + ! call endrun(msg=errMsg(sourcefile, __LINE__)) + ! end if + !end if + !agross = min(r1,r2) - agross = min(r1,r2) + agross = QuadraticSolveSmooth_theta_ip(bquad, cquad) end if @@ -1446,17 +1539,17 @@ function LeafHumidityStomaResis(leaf_psi, k_lwp, veg_tempk, can_vpress, can_pres ! ------------------------------------------------------------------------------------- ! Arguments - real(r8) :: leaf_psi ! Leaf water potential [MPa] - real(r8) :: k_lwp ! Scaling coefficient for the ratio of leaf xylem (user parameter) - real(r8) :: veg_tempk ! Leaf temperature [K] - real(r8) :: can_vpress ! vapor pressure of air (unconstrained) [Pa] - real(r8) :: can_press ! Atmospheric pressure of canopy [Pa] - - real(r8) :: rb ! Leaf Boundary layer resistance [s/m] - real(r8) :: gstoma ! Stomatal Conductance of this leaf layer [m/s] - integer :: ft ! Plant Functional Type + real(r8),intent(in) :: leaf_psi ! Leaf water potential [MPa] + real(r8),intent(in) :: k_lwp ! Scaling coefficient for the ratio of leaf xylem (user parameter) + real(r8),intent(in) :: veg_tempk ! Leaf temperature [K] + real(r8),intent(in) :: can_vpress ! vapor pressure of air (unconstrained) [Pa] + real(r8),intent(in) :: can_press ! Atmospheric pressure of canopy [Pa] + real(r8),intent(in) :: rb ! Leaf Boundary layer resistance [s/m] + real(r8),intent(in) :: gstoma ! Stomatal Conductance of this leaf layer [m/s] + integer,intent(in) :: ft ! Plant Functional Type + real(r8),intent(in) :: veg_esat ! Saturation vapor pressure at veg surface [Pa] + real(r8) :: rstoma_out ! Total Stomatal resistance (stoma and BL) [s/m] - real(r8) :: veg_esat ! Saturation vapor pressure at veg surface [Pa] ! Locals real(r8) :: ceair ! vapor pressure of air, constrained [Pa] @@ -1545,7 +1638,6 @@ function ft1_f(tl, ha) result(ans) ! 7/23/16: Copied over from CLM by Ryan Knox ! !!USES - ! ! !ARGUMENTS: @@ -1562,7 +1654,41 @@ function ft1_f(tl, ha) result(ans) end function ft1_f ! ===================================================================================== + + elemental function ft1_fo(tl, ft1_a, ft1_b) result(ans) + ! + !!DESCRIPTION: + ! optimized photosynthesis temperature response + ! + ! This function, along with the fth1_f() function, scales the Vcmax25 and Jmax25 + ! constants so they are adjusted for leaf temperature. + ! This equation does not have a name, but is equation 9.10 in the CLM5.0 technical + ! manual. + ! + ! !REVISION HISTORY + ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 + ! 7/23/16: Copied over from CLM by Ryan Knox + ! 3/26: Refactored for speed + ! + !!USES + ! + ! !ARGUMENTS: + real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temperature function (K) + real(r8), intent(in) :: ft1_a + real(r8), intent(in) :: ft1_b + ! + ! !LOCAL VARIABLES: + real(r8) :: ans + !------------------------------------------------------------------------------- + + ans = exp( ft1_a - ft1_b/tl) + + return + end function ft1_fo + + ! ===================================================================================== + function fth_f(tl,hd,se,scaleFactor) result(ans) ! !!DESCRIPTION: @@ -1597,6 +1723,38 @@ end function fth_f ! ===================================================================================== + elemental function fth_fo(tl,fth_a,fth_b,scale_factor) result(ans) + ! + !!DESCRIPTION: + ! optimized photosynthesis temperature inhibition + ! + ! This function, along with the ft1_f() function, scales the Vcmax25 and Jmax25 + ! constants so they are adjusted for leaf temperature. + ! This equation does not have a name, but is equation 9.10 in the CLM5.0 technical + ! manual. + ! !REVISION HISTORY + ! Jinyun Tang separated it out from Photosynthesis, Feb. 07/2013 + ! 7/23/16: Copied over from CLM by Ryan Knox + ! + + ! + ! !ARGUMENTS: + real(r8), intent(in) :: tl ! leaf temperature in photosynthesis temp function (K) + real(r8), intent(in) :: fth_a + real(r8), intent(in) :: fth_b + real(r8), intent(in) :: scale_factor ! scaling factor for high temp inhibition (25 C = 1.0) + ! + ! !LOCAL VARIABLES: + real(r8) :: ans + !------------------------------------------------------------------------------- + + ans = scale_factor / ( 1._r8 + exp( fth_a + fth_b/tl ) ) + + return + end function fth_fo + + ! ===================================================================================== + function fth25_f(hd,se)result(ans) ! !!DESCRIPTION: @@ -1708,11 +1866,9 @@ subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, & nscaler, & ft, & veg_tempk, & + bprate, & lmr) - - - ! ----------------------------------------------------------------------- ! Base maintenance respiration rate for plant tissues maintresp_leaf_ryan1991_baserate ! M. Ryan, 1991. Effects of climate change on plant respiration. @@ -1726,39 +1882,24 @@ subroutine LeafLayerMaintenanceRespiration_Ryan_1991(lnc_top, & real(r8), intent(in) :: nscaler ! Scale for leaf nitrogen profile integer, intent(in) :: ft ! (plant) Functional Type Index real(r8), intent(in) :: veg_tempk ! vegetation temperature + type(bprate_type),intent(in) :: bprate! Temperature response function 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) :: lmr25top ! canopy top leaf maint resp rate at 25C for this pft (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) - - lmr25top = lb_params%maintresp_leaf_ryan1991_baserate(ft) * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) - lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) - + + ! Baserate is in umolC/gN/s @25C + ! umolC/gN/s*gN/m2 = umolC/m2/s + lmr25top = lb_params%maintresp_leaf_ryan1991_baserate(ft) * lnc_top ! Part I: Leaf Maintenance respiration: umol CO2 / m**2 [leaf] / s ! ---------------------------------------------------------------------------------- lmr25 = lmr25top * nscaler - - if (lb_params%c3psn(ft) == c3_path_index) then - ! temperature sensitivity of C3 plants - lmr = lmr25 * ft1_f(veg_tempk, lmrha) * & - fth_f(veg_tempk, lmrhd, lmrse, lmrc) - else - ! temperature sensitivity of C4 plants - lmr = lmr25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) - lmr = lmr / (1._r8 + exp( 1.3_r8*(veg_tempk-(tfrz+55._r8)) )) - endif - + ! temperature response + lmr = lmr25* bprate%lmr_tscaler(ft) + ! Any hydrodynamic limitations could go here, currently none ! lmr = lmr * (nothing) @@ -1771,16 +1912,16 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & ft, & veg_tempk, & tgrowth, & + bprate, & lmr) - - ! Arguments real(r8), intent(in) :: lnc_top ! Leaf nitrogen content per unit area at canopy top [gN/m2] real(r8), intent(in) :: rdark_scaler ! Decay coefficient for vertical scaling of respiration. See note 1 below integer, intent(in) :: ft ! (plant) Functional Type Index real(r8), intent(in) :: veg_tempk ! vegetation temperature (degrees K) real(r8), intent(in) :: tgrowth ! lagged vegetation temperature averaged over acclimation timescale (degrees K) + type(bprate_type),intent(in) :: bprate ! Temperature response function real(r8), intent(out) :: lmr ! Leaf Maintenance Respiration (umol CO2/m**2/s) @@ -1816,8 +1957,10 @@ subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017(lnc_top, & call FatesWarn(warn_msg,index=4) end if - lmr = r_t_ref * exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * & - ((veg_tempk-tfrz)**2 - lmr_TrefC**2)) + lmr = r_t_ref * bprate%lmr_tscaler(ft) + + !exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * & + ! ((veg_tempk-tfrz)**2._r8 - lmr_TrefC**2._r8)) end subroutine LeafLayerMaintenanceRespiration_Atkin_etal_2017 @@ -1830,8 +1973,7 @@ subroutine LeafLayerBiophysicalRates( ft, & nscaler, & veg_tempk, & dayl_factor, & - t_growth, & - t_home, & + bprate, & btran, & vcmax, & jmax, & @@ -1857,28 +1999,27 @@ subroutine LeafLayerBiophysicalRates( ft, & ! Arguments ! ------------------------------------------------------------------------------ - 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) :: kp25_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) :: dayl_factor ! daylength scaling factor (0-1) - real(r8), intent(in) :: t_growth ! T_growth (short-term running mean temperature) (K) - real(r8), intent(in) :: t_home ! T_home (long-term running mean temperature) (K) - 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) :: kp ! initial slope of CO2 response curve (C4 plants) - real(r8), intent(out) :: gs0 ! effective stomatal intercept - real(r8), intent(out) :: gs1 ! effective stomatal slope - real(r8), intent(out) :: gs2 ! alternative btran term applied to Medlyn - ! conductance on whole non-intercept side of equation - + 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) :: kp25_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) :: dayl_factor ! daylength scaling factor (0-1) + type(bprate_type),intent(in) :: bprate + 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) :: kp ! initial slope of CO2 response curve (C4 plants) + real(r8), intent(out) :: gs0 ! effective stomatal intercept + real(r8), intent(out) :: gs1 ! effective stomatal slope + real(r8), intent(out) :: gs2 ! alternative btran term applied to Medlyn + ! conductance on whole non-intercept side of equation + ! Locals ! ------------------------------------------------------------------------------- real(r8) :: vcmax25 ! leaf layer: maximum rate of carboxylation at 25C @@ -1889,18 +2030,6 @@ subroutine LeafLayerBiophysicalRates( ft, & ! Parameters ! --------------------------------------------------------------------------------- - real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) - real(r8) :: jmaxha ! activation energy for jmax (J/mol) - real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) - real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) - real(r8) :: vcmaxse ! entropy term for vcmax (J/mol/K) - real(r8) :: jmaxse ! entropy term for jmax (J/mol/K) - real(r8) :: t_growth_celsius ! average growing temperature - real(r8) :: t_home_celsius ! average home temperature - real(r8) :: jvr ! ratio of Jmax25 / Vcmax25 - 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) - ! update the daylength factor local variable if the switch is on if ( lb_params%dayl_switch == itrue ) then @@ -1909,36 +2038,6 @@ subroutine LeafLayerBiophysicalRates( ft, & dayl_factor_local = 1.0_r8 endif - select case(lb_params%photo_tempsens_model) - case (photosynth_acclim_model_none) !No temperature acclimation - vcmaxha = lb_params%vcmaxha(FT) - jmaxha = lb_params%jmaxha(FT) - vcmaxhd = lb_params%vcmaxhd(FT) - jmaxhd = lb_params%jmaxhd(FT) - vcmaxse = lb_params%vcmaxse(FT) - jmaxse = lb_params%jmaxse(FT) - - case (photosynth_acclim_model_kumarathunge_etal_2019) !Kumarathunge et al. temperature acclimation, Thome=30-year running mean - t_growth_celsius = t_growth-tfrz - t_home_celsius = t_home-tfrz - vcmaxha = (42.6_r8 + (1.14_r8*t_growth_celsius))*1e3_r8 !J/mol - jmaxha = 40.71_r8*1e3_r8 !J/mol - vcmaxhd = 200._r8*1e3_r8 !J/mol - jmaxhd = 200._r8*1e3_r8 !J/mol - vcmaxse = (645.13_r8 - (0.38_r8*t_growth_celsius)) - jmaxse = 658.77_r8 - (0.84_r8*t_home_celsius) - 0.52_r8*(t_growth_celsius-t_home_celsius) - jvr = 2.56_r8 - (0.0375_r8*t_home_celsius)-(0.0202_r8*(t_growth_celsius-t_home_celsius)) - - - case default - write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' - write (fates_log(),*)'lb_params%photo_tempsens_model: ',lb_params%photo_tempsens_model - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select - - vcmaxc = fth25_f(vcmaxhd, vcmaxse) - jmaxc = fth25_f(jmaxhd, jmaxse) - ! Vcmax25top was already calculated to derive the nscaler function vcmax25 = vcmax25top_ft * nscaler * dayl_factor_local @@ -1946,7 +2045,7 @@ subroutine LeafLayerBiophysicalRates( ft, & case (photosynth_acclim_model_none) jmax25 = jmax25top_ft * nscaler * dayl_factor_local case (photosynth_acclim_model_kumarathunge_etal_2019) - jmax25 = vcmax25*jvr + jmax25 = vcmax25*bprate%jvr(ft) case default write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' write (fates_log(),*)'lb_params%photo_tempsens_model:',lb_params%photo_tempsens_model @@ -1957,7 +2056,7 @@ subroutine LeafLayerBiophysicalRates( ft, & ! photosynthetic pathway: 0. = c4, 1. = c3 if (lb_params%c3psn(ft) == c3_path_index) then - vcmax = vcmax25 * ft1_f(veg_tempk, vcmaxha) * fth_f(veg_tempk, vcmaxhd, vcmaxse, vcmaxc) + vcmax = vcmax25 * bprate%vcmax_tscaler(ft) kp = -9999._r8 else vcmax = vcmax25 * 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) @@ -1966,8 +2065,8 @@ subroutine LeafLayerBiophysicalRates( ft, & kp = kp25_ft * nscaler * 2._r8**((min(veg_tempk,310._r8)-(tfrz+25._r8))/10._r8) end if - jmax = jmax25 * ft1_f(veg_tempk, jmaxha) * fth_f(veg_tempk, jmaxhd, jmaxse, jmaxc) - + jmax = jmax25 * bprate%jmax_tscaler(ft) + ! Adjust various rates for water limitations ! ----------------------------------------------------------------------------------- @@ -2028,9 +2127,6 @@ subroutine LeafLayerBiophysicalRates( ft, & gs1 = lb_params%bb_slope(ft) end if end if - - - return end subroutine LeafLayerBiophysicalRates @@ -2121,7 +2217,7 @@ end subroutine LowstorageMainRespReduction ! ===================================================================================== - real(r8) function GetConstrainedVPress(air_vpress,veg_esat) result(ceair) + function GetConstrainedVPress(air_vpress,veg_esat) result(ceair) ! ----------------------------------------------------------------------------------- ! Return a constrained vapor pressure [Pa] @@ -2133,18 +2229,20 @@ real(r8) function GetConstrainedVPress(air_vpress,veg_esat) result(ceair) ! behaved. ! ----------------------------------------------------------------------------------- - real(r8) :: air_vpress ! vapor pressure of the air (unconstrained) [Pa] - real(r8) :: veg_esat ! saturated vapor pressure [Pa] - real(r8) :: min_frac_esat = 0.05_r8 ! We don't allow vapor pressures - ! below this fraction amount of saturation vapor pressure + real(r8),intent(in) :: air_vpress ! vapor pressure of the air (unconstrained) [Pa] + real(r8),intent(in) :: veg_esat ! saturated vapor pressure [Pa] + real(r8) :: ceair - ceair = min( max(air_vpress, min_frac_esat*veg_esat ),veg_esat ) + real(r8),parameter :: min_frac_esat = 0.05_r8 ! We don't allow vapor pressures + ! below this fraction amount of saturation vapor pressure + + ceair = min( max(air_vpress, min_frac_esat*veg_esat),veg_esat ) end function GetConstrainedVPress ! ===================================================================================== - subroutine QSat (T, p, qs, es, qsdT, esdT) + elemental subroutine QSat (T, p, qs, es, qsdT, esdT) ! ! !DESCRIPTION: ! @@ -2257,7 +2355,7 @@ end subroutine QSat ! ===================================================================================== - real(r8) function VeloToMolarCF(press,tempk) result(cf) + function VeloToMolarCF(press,tempk) result(cf) ! --------------------------------------------------------------------------------- ! @@ -2287,8 +2385,9 @@ real(r8) function VeloToMolarCF(press,tempk) result(cf) ! -------------------------------------------------------------------------------- ! Arguments - real(r8) :: press ! air pressure at point of interest [Pa] - real(r8) :: tempk ! temperature at point of interest [K] + real(r8),intent(in) :: press ! air pressure at point of interest [Pa] + real(r8),intent(in) :: tempk ! temperature at point of interest [K] + real(r8) :: cf cf = press/(rgas_J_K_kmol * tempk )*umol_per_kmol @@ -2296,4 +2395,271 @@ end function VeloToMolarCF ! ===================================================================================== + elemental function QuadraticSolveSmooth_a1_max(b, c) result(root) + + ! ----------------------------------------------------------------------------------- + ! Specialized solver for a = 1.0 - returns LARGER root (max) + ! Used for: Medlyn stomatal conductance + ! Returns: (-b + sqrt(b^2 - 4*1*c)) / 2*1 + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: b, c + real(r8) :: root + real(r8) :: discriminant + real(r8), parameter :: inv_2 = 0.5_r8 + + discriminant = b * b - 4.0_r8 * c + + root = (-b + sqrt(max(discriminant, 0.0_r8))) * inv_2 + + end function QuadraticSolveSmooth_a1_max + + ! ===================================================================================== + + elemental function QuadraticSolveSmooth_theta_psii(b, c) result(root) + + ! ----------------------------------------------------------------------------------- + ! Specialized solver for a = theta_psii = 0.7 (electron transport) + ! Returns the SMALLER root (min): (-b - sqrt(b^2 - 4*theta*c)) / (2*theta) + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: b, c + real(r8) :: root + real(r8) :: r1,r2 + real(r8) :: discriminant,dsq + real(r8), parameter :: four_theta_psii = 4.0_r8 * theta_psii + real(r8), parameter :: inv_2_theta_psii = 1.0_r8 / (2.0_r8 * theta_psii) + + discriminant = b * b - four_theta_psii * c + dsq = sqrt(max(discriminant, 0.0_r8)) + !dsq = sqrt(abs(discriminant)) + r1 = (-b - dsq) * inv_2_theta_psii + r2 = (-b + dsq) * inv_2_theta_psii + root = min(r1,r2) + + end function QuadraticSolveSmooth_theta_psii + + ! ===================================================================================== + + elemental function QuadraticSolveSmooth_theta_cj(b, c) result(root) + + ! ----------------------------------------------------------------------------------- + ! Specialized solver for a = theta_cj_c4 = 0.98 (C4 RuBP limitation) + ! Returns the SMALLER root (min): + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: b, c + real(r8) :: root + real(r8) :: discriminant + real(r8), parameter :: four_theta_cj = 4.0_r8 * theta_cj_c4 + real(r8), parameter :: inv_2_theta_cj = 1.0_r8 / (2.0_r8 * theta_cj_c4) + + discriminant = b * b - four_theta_cj * c + + root = (-b - sqrt(max(discriminant, 0.0_r8))) * inv_2_theta_cj + + end function QuadraticSolveSmooth_theta_cj + + ! ===================================================================================== + + elemental function QuadraticSolveSmooth_theta_ip(b, c) result(root) + + ! ----------------------------------------------------------------------------------- + ! Specialized solver for a = theta_ip_c4 = 0.95 (C4 PEP limitation) + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: b, c + real(r8) :: root + real(r8) :: discriminant + real(r8), parameter :: four_theta_ip = 4.0_r8 * theta_ip_c4 + real(r8), parameter :: inv_2_theta_ip = 1.0_r8 / (2.0_r8 * theta_ip_c4) + + discriminant = b * b - four_theta_ip * c + + root = (-b - sqrt(max(discriminant, 0.0_r8))) * inv_2_theta_ip + + end function QuadraticSolveSmooth_theta_ip + + + elemental subroutine QuadraticSolveSmooth(a, b, c, r1, r2) + + ! ----------------------------------------------------------------------------------- + ! Generic fast solver + ! Used for BB, where a is positive. + ! BB will use the larger root. + ! ----------------------------------------------------------------------------------- + + real(r8), intent(in) :: a, b, c + real(r8), intent(out) :: r1,r2 + real(r8) :: discriminant + real(r8) :: dsq ! sqrt(disc) + + if (abs(a) < nearzero ) then + r2 = 0.0_r8 + if ( abs(b)>nearzero ) r2 = -c/b + r1 = 0.0_r8 + return + end if + + discriminant = b * b - 4._r8*a*c + dsq = sqrt(max(discriminant, 0.0_r8)) + + r1 = 0.5_r8*(-b + dsq)/a + r2 = 0.5_r8*(-b - dsq)/a + + end subroutine QuadraticSolveSmooth + + ! ===================================================================================== + + subroutine UpdateFastBiophysicalRates(bprate,ft,veg_tempk,c3_psn,maintresp_leaf_model) + + type(bprate_type),intent(inout) :: bprate + integer ,intent(in) :: ft + real(r8) ,intent(in) :: veg_tempk + integer ,intent(in) :: c3_psn ! 1=C3, 0=C4 + integer ,intent(in) :: maintresp_leaf_model ! Ryan91 or Atkin2017 + + ! 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 + + ! Update the temperature scalers that apply to vcmax25 and jmax25 + + bprate%vcmax_tscaler(ft) = ft1_fo(veg_tempk, bprate%vr_a(ft), bprate%vr_b(ft)) * & + fth_fo(veg_tempk, bprate%vi_a(ft), bprate%vi_b(ft), bprate%vi_c(ft)) + + bprate%jmax_tscaler(ft) = ft1_fo(veg_tempk, bprate%jr_a(ft), bprate%jr_b(ft)) * & + fth_fo(veg_tempk, bprate%ji_a(ft), bprate%ji_b(ft), bprate%ji_c(ft)) + + if(maintresp_leaf_model == lmrmodel_ryan_1991)then + + if (c3_psn == c3_path_index) then + ! LMR temperature sensitivity of C3 plants + bprate%lmr_tscaler(ft) = ft1_f(veg_tempk, lmrha) * & + fth_f(veg_tempk, lmrhd, lmrse, lmrc) + else + ! temperature sensitivity of C4 plants + bprate%lmr_tscaler(ft) = 2._r8**((veg_tempk-(tfrz+25._r8))/10._r8) / & + (1._r8 + exp( 1.3_r8*(veg_tempk-(tfrz+55._r8)) )) + endif + else + bprate%lmr_tscaler(ft) = exp(lmr_b * (veg_tempk - tfrz - lmr_TrefC) + lmr_c * & + ((veg_tempk-tfrz)**2._r8 - lmr_TrefC**2._r8)) + end if + + + return + end subroutine UpdateFastBiophysicalRates + + + ! ===================================================================================== + + subroutine UpdateSlowBiophysicalRates(bprate,ft,t_growth_k,t_home_k) + + use FatesConstantsMod, only: tfrz => t_water_freeze_k_1atm + + type(bprate_type),intent(inout) :: bprate + integer ,intent(in) :: ft + real(r8) ,intent(in) :: t_growth_k + real(r8) ,intent(in) :: t_home_k + real(r8) :: vcmaxha ! activation energy for vcmax (J/mol) + real(r8) :: jmaxha ! activation energy for jmax (J/mol) + real(r8) :: vcmaxhd ! deactivation energy for vcmax (J/mol) + real(r8) :: jmaxhd ! deactivation energy for jmax (J/mol) + real(r8) :: vcmaxse ! entropy term for vcmax (J/mol/K) + real(r8) :: jmaxse ! entropy term for jmax (J/mol/K) + real(r8) :: jvr ! ratio of Jmax25 / Vcmax25 ! Kumarathunge et al. + 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) :: t_growth_celsius + real(r8) :: t_home_celsius + + select case(lb_params%photo_tempsens_model) + case (photosynth_acclim_model_none) !No temperature acclimation + vcmaxha = lb_params%vcmaxha(ft) + jmaxha = lb_params%jmaxha(ft) + vcmaxhd = lb_params%vcmaxhd(ft) + jmaxhd = lb_params%jmaxhd(ft) + vcmaxse = lb_params%vcmaxse(ft) + jmaxse = lb_params%jmaxse(ft) + case (photosynth_acclim_model_kumarathunge_etal_2019) !Kumarathunge et al. temperature acclimation, Thome=30-year running mean + t_growth_celsius = t_growth_k-tfrz + t_home_celsius = t_home_k-tfrz + vcmaxha = (42.6_r8 + (1.14_r8*t_growth_celsius))*1e3_r8 !J/mol + jmaxha = 40.71_r8*1e3_r8 !J/mol + vcmaxhd = 200._r8*1e3_r8 !J/mol + jmaxhd = 200._r8*1e3_r8 !J/mol + vcmaxse = (645.13_r8 - (0.38_r8*t_growth_celsius)) + jmaxse = 658.77_r8 - (0.84_r8*t_home_celsius) - 0.52_r8*(t_growth_celsius-t_home_celsius) + bprate%jvr(ft) = 2.56_r8 - (0.0375_r8*t_home_celsius)-(0.0202_r8*(t_growth_celsius-t_home_celsius)) + case default + write (fates_log(),*)'error, incorrect leaf photosynthesis temperature acclimation model specified' + write (fates_log(),*)'lb_params%photo_tempsens_model: ',lb_params%photo_tempsens_model + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + + ! Response function terms + bprate%vr_a(ft) = vcmaxha / (rgas_J_K_mol*t_ref_298K) + bprate%vr_b(ft) = vcmaxha / rgas_J_K_mol + bprate%jr_a(ft) = jmaxha / (rgas_J_K_mol*t_ref_298K) + bprate%jr_b(ft) = jmaxha / rgas_J_K_mol + + ! Inhibition function terms + bprate%vi_a(ft) = vcmaxse / rgas_J_K_mol + bprate%vi_b(ft) = -vcmaxhd / rgas_J_K_mol + bprate%vi_c(ft) = fth25_f(vcmaxhd, vcmaxse) + + bprate%ji_a(ft) = jmaxse / rgas_J_K_mol + bprate%ji_b(ft) = -jmaxhd / rgas_J_K_mol + bprate%ji_c(ft) = fth25_f(jmaxhd, jmaxse) + + + return + end subroutine UpdateSlowBiophysicalRates + + + + + function FastQ10(temperature,offset,q10_base,log_q10_div10) result(t_factor) + + ! ---------------------------------------------------------------------------------- + ! The exp() intrinsic is faster than a pow() + ! or ** math operation. So we re-write the + ! q10 function leverage this. Differences + ! are on the order of e-15 on gnu compiler + ! + ! identity: + ! since: a^b = e^{b ln(a)} + ! + ! a^(x-c)/b = e^{(x-c)/b ln(a)} + ! = e^{(x-c) ln(a)/b} + ! + ! original: + ! t_factor_orig = q10_mr**((temp - tfrz - 20.0_r8)/10.0_r8) + ! + ! log_q10_mr_div10 = log(q10_mr)/10. + ! ---------------------------------------------------------------------------------- + + real(r8),intent(in) :: temperature ! The temperature dictating the response [K] + real(r8),intent(in) :: offset ! Offset of the temperature response [C] (e.g. 20C,25C,etc) + real(r8),intent(in) :: q10_base + real(r8),intent(in) :: log_q10_div10 ! = ln(a)/10 + real(r8) :: t_factor + + if(do_b4b)then + t_factor = q10_base**((temperature - tfrz - offset)/10.0_r8) + + else + t_factor = exp(log_q10_div10 * (temperature - tfrz - offset)) + end if + + return + end function FastQ10 + + + + + end module LeafBiophysicsMod diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 684ae8fb6c..0b11877006 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -67,6 +67,7 @@ module EDInitMod use FatesInterfaceTypesMod , only : nlevdamage use FatesInterfaceTypesMod , only : hlm_use_nocomp use FatesInterfaceTypesMod , only : nlevage + use EDParamsMod , only : maxpatch_total use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : bagw_allom @@ -164,6 +165,9 @@ subroutine init_site_vars( site_in, bc_in, bc_out ) allocate(site_in%mass_balance(1:num_elements)) allocate(site_in%iflux_balance(1:num_elements)) + ! Patch type vector + allocate(site_in%pa_vec(maxpatch_total)) + if (hlm_use_tree_damage .eq. itrue) then allocate(site_in%term_nindivs_canopy_damage(1:nlevdamage, 1:nlevsclass, 1:numpft)) allocate(site_in%term_nindivs_ustory_damage(1:nlevdamage, 1:nlevsclass, 1:numpft)) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index f07a3c75de..42ae84848d 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -879,11 +879,13 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) ! Update the total area of by patch age class array currentSite%area_by_age(currentPatch%age_class) = & currentSite%area_by_age(currentPatch%age_class) + currentPatch%area - + currentPatch => currentPatch%younger - enddo + ! Only need to do this if parallel... + call PatchLoads(currentSite,bc_out) + ! Check to see if the time integrated fluxes match the state ! Dont call this if we are restarting, it will double count the flux if(.not.is_restarting)then @@ -921,8 +923,98 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) bc_out%fire_closs_to_atm_si = sum(site_cmass%burn_flux_to_atm(:)) * area_inv * days_per_sec bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * area_inv * days_per_sec + + + end subroutine ed_update_site + + subroutine PatchLoads(site,bc_out) + + + use FatesRadiationMemMod, only: norman_solver + use FatesInterfaceTypesMod, only: hlm_radiation_model + + ! ----------------------------------------------------------------- + ! We calculate the computational burden of a patch to + ! help with multithreading and balancing the load of + ! of threads, which is organized at the patch level + ! Multithreading is applied to photosynthesis and radiation + ! currently. + ! + ! Total photosynthesis calculations (without hydro) are proportional + ! to the total number of leaf layers unique to canopy layer and pft + ! + ! Total photosynthesis calculations (with hydro) are proportional + ! to the total number of leaf layers on every single cohort. + ! + ! Respiration calculations and some accounting is also applied at + ! the cohort level + ! + ! Radiation scattering calculations are proportional to the total + ! number of cohorts^3 (currently) + ! + ! Photosynthesis takes about 3x more computation than radiation. + ! ------------------------------------------------------------------ + + type(ed_site_type), intent(in) :: site + type(bc_out_type), intent(inout) :: bc_out + + integer :: npatches + type(fates_patch_type), pointer :: patch + integer :: ifp + + real(r8), parameter :: radfrac = 0.25_r8 + + npatches = site%youngest_patch%patchno + bc_out%load_size(:) = 0._r8 + block + + real(r8) :: psn_load(npatches) + real(r8) :: rad_load(npatches) + + psn_load(:) = 0._r8 + rad_load(:) = 0._r8 + + patch => site%oldest_patch + do while(associated(patch)) + if(patch%nocomp_pft_label .ne. nocomp_bareground)then + ifp = patch%patchno + + if(hlm_use_planthydro.eq.itrue)then + psn_load(ifp) = patch%num_cohorts + else + psn_load(ifp) = sum(patch%nleaf) + end if + + if(hlm_radiation_model.eq.norman_solver)then + rad_load(ifp) = sum(patch%nleaf)**2._r8 + else + rad_load(ifp) = patch%num_cohorts**2._r8 + end if + + end if + patch=>patch%younger + end do + + ! Normalize + if(sum(psn_load)>0._r8)then + psn_load = psn_load/sum(psn_load) + rad_load = rad_load/sum(rad_load) + end if + + patch => site%oldest_patch + do while(associated(patch)) + if(patch%nocomp_pft_label .ne. nocomp_bareground)then + ifp = patch%patchno + bc_out%load_size(ifp) = (1._r8 - radfrac)*psn_load(ifp) + radfrac*rad_load(ifp) + end if + patch=>patch%younger + end do + end block + + end subroutine PatchLoads + !-------------------------------------------------------------------------------! subroutine TotalBalanceCheck (currentSite, call_index, is_restarting ) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 369eb6b9ec..b6b1819d6d 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -7,6 +7,7 @@ module EDParamsMod use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : nearzero use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun use FatesConstantsMod, only : fates_unset_r8 @@ -93,6 +94,9 @@ module EDParamsMod real(r8),protected,public :: q10_mr ! Q10 for respiration rate (for soil fragmenation and plant respiration) (unitless) real(r8),protected,public :: q10_froz ! Q10 for frozen-soil respiration rates (for soil fragmentation) (unitless) + real(r8),protected,public :: log_q10_mr_div10 != log(q10_mr)/10.0_r8 + real(r8),protected,public :: log_q10_froz_div10 ! + ! grazing parameters real(r8),protected,public :: landuse_grazing_carbon_use_eff real(r8),protected,public :: landuse_grazing_maxheight @@ -191,11 +195,12 @@ module EDParamsMod public :: TransferParamsGeneric public :: FatesReportParams public :: GetNVegLayers - contains + ! ==================================================================================== + function GetNVegLayers(treevai) result(nv) real(r8) :: treevai ! The LAI+SAI of the cohort (m2/m2) @@ -381,6 +386,10 @@ subroutine TransferParamsGeneric(pstruct) param_p => pstruct%GetParamFromName("fates_q10_froz") q10_froz = param_p%r_data_scalar + + ! pre-logged for computational efficiency + log_q10_mr_div10 = log(q10_mr)/10.0_r8 + log_q10_froz_div10 = log(q10_froz)/10.0_r8 param_p => pstruct%GetParamFromName("fates_history_sizeclass_bin_edges") allocate(ED_val_history_sizeclass_bin_edges(size(param_p%r_data_1d,dim=1))) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 5aca294987..4aeba01c37 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -31,6 +31,7 @@ module EDTypesMod use FatesInterfaceTypesMod,only : bc_out_type use FatesConstantsMod , only : n_landuse_cats use FatesInterfaceTypesMod,only : hlm_parteh_mode + use EDParamsMod, only : maxpatch_total use FatesCohortMod, only : fates_cohort_type use FatesPatchMod, only : fates_patch_type use EDParamsMod, only : nclmax, nlevleaf, maxpft @@ -318,15 +319,33 @@ module EDTypesMod end type site_massbal_type + type :: fates_patch_vec_type + + ! This is a scratch array for patch pointers + ! this is useful if you want to loop over patches + ! in order by index. This should be updated + ! upon canopy summarization, and should + ! facilitate coupling with the host + + type(fates_patch_type), pointer :: p => null() + + end type fates_patch_vec_type + + !************************************ !** Site type structure ** !************************************ + + + type, public :: ed_site_type ! POINTERS type (fates_patch_type), pointer :: oldest_patch => null() ! pointer to oldest patch at the site type (fates_patch_type), pointer :: youngest_patch => null() ! pointer to yngest patch at the site + + type (fates_patch_vec_type), allocatable :: pa_vec(:) ! Patch vector by patch no ! Resource management type (ed_resources_management_type) :: resources_management ! resources_management at the site @@ -623,13 +642,18 @@ subroutine set_patchno( currentSite, check , call_id) ! Special case: For no-comp runs, we treat the bare-ground ! patch as index 0. - type(ed_site_type),intent(in) :: currentSite + type(ed_site_type),intent(inout) :: currentSite logical,intent(in) :: check ! If true, we are checking order, not setting integer,intent(in) :: call_id ! An index used for testing type(fates_patch_type), pointer :: currentPatch integer patchno - + !--------------------------------------------------------------------- + + ! Flush the patch pointers + do patchno = 1,maxpatch_total + currentSite%pa_vec(patchno)%p => null() + end do patchno = 1 currentPatch => currentSite%oldest_patch @@ -647,6 +671,7 @@ subroutine set_patchno( currentSite, check , call_id) call endrun(msg=errMsg(sourcefile, __LINE__)) end if currentPatch%patchno = patchno + currentSite%pa_vec(patchno)%p => currentPatch patchno = patchno + 1 endif currentPatch => currentPatch%younger diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 4d535f9de4..d98a33d7d8 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -322,10 +322,13 @@ module FatesConstantsMod ! freezing point of water at 1 atm (K) real(fates_r8), parameter, public :: t_water_freeze_k_1atm = 273.15_fates_r8 - + ! freezing point of water at triple point (K) real(fates_r8), parameter, public :: t_water_freeze_k_triple = 273.16_fates_r8 + ! Reference temperature (at 25C) used in various physiology calculations (K) + real(fates_r8), parameter, public :: t_ref298 = 298.15_fates_r8 + ! Density of fresh liquid water (kg/m3) real(fates_r8), parameter, public :: dens_fresh_liquid_water = 1.0E3_fates_r8 diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 22e207cffc..b189013ffb 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2701,7 +2701,7 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites) hio_elai_si(io_si) = hio_elai_si(io_si) + sum( cpatch%canopy_area_profile(:,:,:) * cpatch%elai_profile(:,:,:) ) * & cpatch%total_canopy_area * AREA_INV - hio_ncl_si(io_si) = hio_ncl_si(io_si) + cpatch%ncl_p * cpatch%area * AREA_INV + hio_ncl_si(io_si) = hio_ncl_si(io_si) + cpatch%ncl * cpatch%area * AREA_INV ! only valid when "strict ppa" enabled if ( comp_excln_exp .lt. 0._r8 ) then @@ -4856,7 +4856,7 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites) end do hio_ncl_si_age(io_si,cpatch%age_class) = hio_ncl_si_age(io_si,cpatch%age_class) & - + cpatch%ncl_p * patch_area_div_site_area + + cpatch%ncl * patch_area_div_site_area hio_fracarea_burnt_si_age(io_si,cpatch%age_class) = hio_fracarea_burnt_si_age(io_si,cpatch%age_class) + & cpatch%frac_burnt / sec_per_day & ! [frac/day] -> [frac/sec] @@ -5172,6 +5172,7 @@ subroutine update_history_hifrq_sitelevel(this,nc,nsites,sites,bc_in,dt_tstep) ! Locals integer :: s ! The local site index integer :: io_si ! The site index of the IO array + integer :: ico ! cohort index integer :: age_class ! class age index real(r8) :: site_area_veg_inv ! inverse canopy area of the site (1/m2) real(r8) :: site_area_rad_inv ! inverse canopy area of site for only @@ -5314,9 +5315,10 @@ subroutine update_history_hifrq_sitelevel(this,nc,nsites,sites,bc_in,dt_tstep) cpatch%total_canopy_area * site_area_veg_inv end if - ccohort => cpatch%shortest + ico = 0 + ccohort => cpatch%tallest do while(associated(ccohort)) - + ico = ico+1 n_perm2 = ccohort%n * AREA_INV if_notnew: if ( .not. ccohort%isnew ) then @@ -5326,44 +5328,44 @@ subroutine update_history_hifrq_sitelevel(this,nc,nsites,sites,bc_in,dt_tstep) ! Net Ecosystem Production [kgC/m2/s]. Use yesterday's growth respiration hio_nep_si(io_si) = hio_nep_si(io_si) + & - (ccohort%gpp_tstep-ccohort%resp_m_tstep) * n_perm2 * dt_tstep_inv - & + (cpatch%coarrays%gpp_tstep(ico)-cpatch%coarrays%resp_m_tstep(ico)) * n_perm2 * dt_tstep_inv - & (ccohort%resp_g_acc_hold+ccohort%resp_excess_hold) * n_perm2 / days_per_year / sec_per_day hio_gpp_si(io_si) = hio_gpp_si(io_si) + & - ccohort%gpp_tstep * n_perm2 * dt_tstep_inv + cpatch%coarrays%gpp_tstep(ico) * n_perm2 * dt_tstep_inv hio_maint_resp_si(io_si) = hio_maint_resp_si(io_si) + & - ccohort%resp_m_tstep * n_perm2 * dt_tstep_inv - + cpatch%coarrays%resp_m_tstep(ico) * n_perm2 * dt_tstep_inv + hio_maint_resp_unreduced_si(io_si) = hio_maint_resp_unreduced_si(io_si) + & - ccohort%resp_m_unreduced * n_perm2 * dt_tstep_inv + cpatch%coarrays%resp_m_unreduced(ico) * n_perm2 * dt_tstep_inv ! Maintenance respiration of different organs - hio_leaf_mr_si(io_si) = hio_leaf_mr_si(io_si) + ccohort%rdark & - * n_perm2 - hio_froot_mr_si(io_si) = hio_froot_mr_si(io_si) + ccohort%froot_mr & - * n_perm2 - hio_livecroot_mr_si(io_si) = hio_livecroot_mr_si(io_si) + ccohort%livecroot_mr & - * n_perm2 - hio_livestem_mr_si(io_si) = hio_livestem_mr_si(io_si) + ccohort%livestem_mr & - * n_perm2 - + hio_leaf_mr_si(io_si) = hio_leaf_mr_si(io_si) + & + cpatch%coarrays%rdark_tstep(ico) * n_perm2 + hio_froot_mr_si(io_si) = hio_froot_mr_si(io_si) + & + cpatch%coarrays%froot_mr(ico) * n_perm2 + hio_livecroot_mr_si(io_si) = hio_livecroot_mr_si(io_si) + & + cpatch%coarrays%livecroot_mr(ico) * n_perm2 + hio_livestem_mr_si(io_si) = hio_livestem_mr_si(io_si) + & + cpatch%coarrays%livestem_mr(ico) * n_perm2 + ! accumulate fluxes on canopy- and understory- separated fluxes ! these fluxes have conversions of [kg/plant/timestep] -> [kg/m2/s] if (ccohort%canopy_layer .eq. 1) then hio_gpp_canopy_si(io_si) = hio_gpp_canopy_si(io_si) + & - ccohort%gpp_tstep * n_perm2 * dt_tstep_inv + cpatch%coarrays%gpp_tstep(ico) * n_perm2 * dt_tstep_inv else hio_gpp_understory_si(io_si) = hio_gpp_understory_si(io_si) + & - ccohort%gpp_tstep * n_perm2 * dt_tstep_inv + cpatch%coarrays%gpp_tstep(ico) * n_perm2 * dt_tstep_inv end if end if if_notnew - ccohort => ccohort%taller + ccohort => ccohort%shorter end do cpatch => cpatch%younger end do @@ -5403,6 +5405,7 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) integer :: lb1,ub1,lb2,ub2 ! IO array bounds for the calling thread integer :: ivar ! index of IO variable object vector integer :: ft ! functional type index + integer :: ico real(r8) :: n_density ! individual of cohort per m2. real(r8) :: n_perm2 ! individuals per m2 for the whole column real(r8) :: site_area_veg_inv ! 1/area of the site that is not bare-ground @@ -5485,10 +5488,11 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) patch_loop1: do while(associated(cpatch)) nocomp_bare: if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then - - ccohort => cpatch%shortest - do while(associated(ccohort)) + ico=0 + ccohort => cpatch%tallest + do while(associated(ccohort)) + ico=ico+1 n_perm2 = ccohort%n * AREA_INV if ( .not. ccohort%isnew ) then @@ -5499,7 +5503,7 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) ! Total AR (kgC/m2/s) = (kgC/plant/step) / (s/step) * (plant/m2) hio_ar_si_scpf(io_si,scpf) = hio_ar_si_scpf(io_si,scpf) + & - (ccohort%resp_m_tstep*dt_tstep_inv) * n_perm2 + & + (cpatch%coarrays%resp_m_tstep(ico)*dt_tstep_inv) * n_perm2 + & (ccohort%resp_g_acc_hold + ccohort%resp_excess_hold)* n_perm2 / days_per_year / sec_per_day ! Growth AR (kgC/m2/s) ! CDK: this should be daily @@ -5508,28 +5512,28 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) ! Maint AR (kgC/m2/s) hio_ar_maint_si_scpf(io_si,scpf) = hio_ar_maint_si_scpf(io_si,scpf) + & - (ccohort%resp_m_tstep*dt_tstep_inv) * n_perm2 + (cpatch%coarrays%resp_m_tstep(ico)*dt_tstep_inv) * n_perm2 ! Maintenance AR partition variables are stored as rates (kgC/plant/s) ! (kgC/m2/s) = (kgC/plant/s) * (plant/m2) hio_ar_agsapm_si_scpf(io_si,scpf) = hio_ar_agsapm_si_scpf(io_si,scpf) + & - ccohort%livestem_mr * n_perm2 + cpatch%coarrays%livestem_mr(ico) * n_perm2 ! (kgC/m2/s) = (kgC/plant/s) * (plant/m2) hio_ar_darkm_si_scpf(io_si,scpf) = hio_ar_darkm_si_scpf(io_si,scpf) + & - ccohort%rdark * n_perm2 + cpatch%coarrays%rdark_tstep(ico) * n_perm2 ! (kgC/m2/s) = (kgC/plant/s) * (plant/m2) hio_ar_crootm_si_scpf(io_si,scpf) = hio_ar_crootm_si_scpf(io_si,scpf) + & - ccohort%livecroot_mr * n_perm2 + cpatch%coarrays%livecroot_mr(ico) * n_perm2 ! (kgC/m2/s) = (kgC/plant/s) * (plant/m2) hio_ar_frootm_si_scpf(io_si,scpf) = hio_ar_frootm_si_scpf(io_si,scpf) + & - ccohort%froot_mr * n_perm2 + cpatch%coarrays%froot_mr(ico) * n_perm2 if (cpatch%land_use_label .gt. nocomp_bareground_land) then hio_gpp_si_landuse(io_si,cpatch%land_use_label) = hio_gpp_si_landuse(io_si,cpatch%land_use_label) & - + ccohort%gpp_tstep * ccohort%n * dt_tstep_inv + + cpatch%coarrays%gpp_tstep(ico) * ccohort%n * dt_tstep_inv end if ! accumulate fluxes on canopy- and understory- separated fluxes @@ -5537,32 +5541,32 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) ! size-resolved respiration fluxes are in kg C / m2 / s hio_rdark_canopy_si_scls(io_si,scls) = hio_rdark_canopy_si_scls(io_si,scls) + & - ccohort%rdark * ccohort%n * ha_per_m2 + cpatch%coarrays%rdark_tstep(ico) * ccohort%n * ha_per_m2 hio_livestem_mr_canopy_si_scls(io_si,scls) = hio_livestem_mr_canopy_si_scls(io_si,scls) + & - ccohort%livestem_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%livestem_mr(ico) * ccohort%n * ha_per_m2 hio_livecroot_mr_canopy_si_scls(io_si,scls) = hio_livecroot_mr_canopy_si_scls(io_si,scls) + & - ccohort%livecroot_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%livecroot_mr(ico) * ccohort%n * ha_per_m2 hio_froot_mr_canopy_si_scls(io_si,scls) = hio_froot_mr_canopy_si_scls(io_si,scls) + & - ccohort%froot_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%froot_mr(ico) * ccohort%n * ha_per_m2 hio_resp_g_canopy_si_scls(io_si,scls) = hio_resp_g_canopy_si_scls(io_si,scls) + & ccohort%resp_g_acc_hold * n_perm2 / days_per_year / sec_per_day hio_resp_m_canopy_si_scls(io_si,scls) = hio_resp_m_canopy_si_scls(io_si,scls) + & - ccohort%resp_m_tstep * ccohort%n * dt_tstep_inv * ha_per_m2 + cpatch%coarrays%resp_m_tstep(ico) * ccohort%n * dt_tstep_inv * ha_per_m2 else ! size-resolved respiration fluxes are in kg C / m2 / s hio_rdark_understory_si_scls(io_si,scls) = hio_rdark_understory_si_scls(io_si,scls) + & - ccohort%rdark * ccohort%n * ha_per_m2 + cpatch%coarrays%rdark_tstep(ico) * ccohort%n * ha_per_m2 hio_livestem_mr_understory_si_scls(io_si,scls) = hio_livestem_mr_understory_si_scls(io_si,scls) + & - ccohort%livestem_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%livestem_mr(ico) * ccohort%n * ha_per_m2 hio_livecroot_mr_understory_si_scls(io_si,scls) = hio_livecroot_mr_understory_si_scls(io_si,scls) + & - ccohort%livecroot_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%livecroot_mr(ico) * ccohort%n * ha_per_m2 hio_froot_mr_understory_si_scls(io_si,scls) = hio_froot_mr_understory_si_scls(io_si,scls) + & - ccohort%froot_mr * ccohort%n * ha_per_m2 + cpatch%coarrays%froot_mr(ico) * ccohort%n * ha_per_m2 hio_resp_g_understory_si_scls(io_si,scls) = hio_resp_g_understory_si_scls(io_si,scls) + & ccohort%resp_g_acc_hold * n_perm2 / days_per_year / sec_per_day hio_resp_m_understory_si_scls(io_si,scls) = hio_resp_m_understory_si_scls(io_si,scls) + & - ccohort%resp_m_tstep * ccohort%n * dt_tstep_inv * ha_per_m2 + cpatch%coarrays%resp_m_tstep(ico) * ccohort%n * dt_tstep_inv * ha_per_m2 endif end associate endif @@ -5572,10 +5576,10 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) do ileaf=1,ccohort%nv cnlf_indx = ileaf + (ican-1) * nlevleaf hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) = hio_ts_net_uptake_si_cnlf(io_si, cnlf_indx) + & - ccohort%ts_net_uptake(ileaf) * dt_tstep_inv * ccohort%c_area * area_inv + cpatch%coarrays%ts_net_uptake(ileaf,ico) * dt_tstep_inv * ccohort%c_area * area_inv end do - ccohort => ccohort%taller + ccohort => ccohort%shorter enddo ! cohort loop @@ -5584,7 +5588,7 @@ subroutine update_history_hifrq_subsite(this,nc,nsites,sites,bc_in,dt_tstep) if_zenith1: if( sites(s)%coszen>0._r8 ) then do_pft1: do ipft=1,numpft - do_canlev1: do ican=1,cpatch%ncl_p + do_canlev1: do ican=1,cpatch%ncl do_leaflev1: do ileaf=1,cpatch%nrad(ican,ipft) ! calculate where we are on multiplexed dimensions @@ -5794,7 +5798,7 @@ subroutine update_history_hifrq_subsite_ageclass(this,nsites,sites,dt_tstep) type(fates_cohort_type), pointer :: ccohort type(fates_patch_type), pointer :: cpatch - integer :: s, io_si + integer :: s, io_si, ico real(r8) :: site_canopy_area real(r8) :: dt_tstep_inv ! Time step in frequency units (/s) real(r8) :: patch_canarea_div_site_canarea ! Weighting based on patch canopy area relative to site canopy area @@ -5842,19 +5846,18 @@ subroutine update_history_hifrq_subsite_ageclass(this,nsites,sites,dt_tstep) hio_c_lblayer_si_age(io_si,cpatch%age_class) = 0._r8 end if - ccohort => cpatch%shortest + ico=0 + ccohort => cpatch%tallest do while(associated(ccohort)) - if (ccohort%isnew) then - ccohort => ccohort%taller - cycle - end if - cohort_n_div_site_area = ccohort%n * AREA_INV + ico=ico+1 + if (.not.ccohort%isnew) then + cohort_n_div_site_area = ccohort%n * AREA_INV - hio_gpp_si_age(io_si,cpatch%age_class) = hio_gpp_si_age(io_si,cpatch%age_class) & - + ccohort%gpp_tstep * dt_tstep_inv & + hio_gpp_si_age(io_si,cpatch%age_class) = hio_gpp_si_age(io_si,cpatch%age_class) & + + cpatch%coarrays%gpp_tstep(ico) * dt_tstep_inv & * cohort_n_div_site_area - - ccohort => ccohort%taller + end if + ccohort => ccohort%shorter end do ! cohort loop cpatch => cpatch%younger diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3ca8220ddb..6ab2630f7e 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -385,6 +385,8 @@ subroutine zero_bcs(fates,s) fates%bc_out(s)%rssun_pa(:) = 0.0_r8 fates%bc_out(s)%rssha_pa(:) = 0.0_r8 + + fates%bc_out(s)%load_size(:) = 0._r8 fates%bc_out(s)%albd_parb(:,:) = 1.0_r8 ! zero albedo, soil absorbs all rad fates%bc_out(s)%albi_parb(:,:) = 1.0_r8 ! zero albedo, soil absorbs all rad @@ -633,10 +635,11 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) allocate(bc_out%rssun_pa(maxpatch_total)) allocate(bc_out%rssha_pa(maxpatch_total)) + + ! Load balancing + allocate(bc_out%load_size(maxpatch_total)) ! Canopy Radiation - - allocate(bc_out%albd_parb(maxpatch_total,num_swb)) allocate(bc_out%albi_parb(maxpatch_total,num_swb)) allocate(bc_out%fabd_parb(maxpatch_total,num_swb)) @@ -2405,7 +2408,7 @@ subroutine SeedlingParPatch(cpatch, & par_low_frac = cpatch%total_canopy_area ! Work up through the canopy layers from the bottom layer - do cl = cpatch%NCL_p,max(1,cpatch%NCL_p-1),-1 + do cl = cpatch%ncl,max(1,cpatch%ncl-1),-1 cl_par = 0._r8 cl_area = 0._r8 do ipft = 1,numpft @@ -2431,7 +2434,7 @@ subroutine SeedlingParPatch(cpatch, & ! Since we are working up through the canopy layers from the ground, ! set the par_high to the previous par_low value and update ! the par_low to the new cl_par value - if(cl .lt. cpatch%NCL_p) then + if(cl .lt. cpatch%ncl) then seedling_par_high = seedling_par_low par_high_frac = (1._r8-cl_area) seedling_par_low = cl_par diff --git a/main/FatesInterfaceTypesMod.F90 b/main/FatesInterfaceTypesMod.F90 index c303720792..0eaee28c82 100644 --- a/main/FatesInterfaceTypesMod.F90 +++ b/main/FatesInterfaceTypesMod.F90 @@ -628,6 +628,10 @@ module FatesInterfaceTypesMod ! Shaded canopy resistance [s/m] real(r8), allocatable :: rssha_pa(:) + + ! For load-balancing by patch + real(r8), allocatable :: load_size(:) + ! leaf photosynthesis (umol CO2 /m**2/ s) ! (NOT CURRENTLY USED, PLACE-HOLDER) !real(r8), allocatable :: psncanopy_pa(:) diff --git a/main/FatesParameterDerivedMod.F90 b/main/FatesParameterDerivedMod.F90 index a044c33721..957666fb30 100644 --- a/main/FatesParameterDerivedMod.F90 +++ b/main/FatesParameterDerivedMod.F90 @@ -24,8 +24,6 @@ module FatesParameterDerivedMod 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 @@ -57,9 +55,7 @@ subroutine InitAllocate(this,numpft) integer, intent(in) :: numpft allocate(this%jmax25top(numpft,nleafage)) - allocate(this%tpu25top(numpft,nleafage)) allocate(this%kp25top(numpft,nleafage)) - allocate(this%branch_frac(numpft)) @@ -115,7 +111,6 @@ subroutine Init(this,numpft) ! (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(ft) this%jmax25top(ft,iage) = 1.67_r8 * vcmax25top(ft,iage) - this%tpu25top(ft,iage) = 0.167_r8 * vcmax25top(ft,iage) this%kp25top(ft,iage) = 20000._r8 * vcmax25top(ft,iage) end do diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index d1ce7689eb..3f60929055 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -2814,7 +2814,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ,io_idx_co,cohortsperpatch endif - this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) = cpatch%ncl_p + this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) = cpatch%ncl this%rvars(ir_zstar_pa)%r81d(io_idx_co_1st) = cpatch%zstar if(hlm_use_sp.eq.ifalse)then @@ -3835,7 +3835,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! Set zenith angle info - cpatch%ncl_p = this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) + cpatch%ncl = this%rvars(ir_nclp_pa)%int1d(io_idx_co_1st) cpatch%zstar = this%rvars(ir_zstar_pa)%r81d(io_idx_co_1st) call this%GetRMeanRestartVar(cpatch%tveg24, ir_tveg24_pa, io_idx_co_1st) @@ -4225,9 +4225,6 @@ subroutine update_3dpatch_radiation(this, nsites, sites, bc_out) call twostr%Solve(ib, & ! in normalized_upper_boundary, & ! in 1.0_r8,1.0_r8, & ! in - sites(s)%taulambda_2str, & ! inout (scratch) - sites(s)%omega_2str, & ! inout (scratch) - sites(s)%ipiv_2str, & ! inout (scratch) albd_parb(ifp,ib), & albi_parb(ifp,ib), & currentPatch%rad_error(ib), & diff --git a/main/FatesUtilsMod.F90 b/main/FatesUtilsMod.F90 index bec737b940..4797cd5eb2 100644 --- a/main/FatesUtilsMod.F90 +++ b/main/FatesUtilsMod.F90 @@ -253,7 +253,7 @@ subroutine QuadraticRootsNSWC(a,b,c,root1,root2,err) end subroutine QuadraticRootsNSWC - subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) + pure subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) real(r8),intent(in) :: a , b , c !! coefficients @@ -273,10 +273,8 @@ subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) end if d = b * b - 4._r8 * a * c - das = sqrt(abs(d)) - if (d > nearzero) then - + das = sqrt(d) root1 = (-b + das) / (2._r8 * a) root2 = (-b - das) / (2._r8 * a) @@ -286,11 +284,11 @@ subroutine QuadraticRootsSridharachary(a,b,c,root1,root2,err) root2 = root1 else - write (fates_log(),*)'error, imaginary roots detected in quadratic solve' + !!write (fates_log(),*)'error, imaginary roots detected in quadratic solve' err = .true. ! Disable this endrun and use the return err to track down ! error provenance - call endrun(msg=errMsg(sourcefile, __LINE__)) + !!call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/parteh/PRTAllometricCNPMod.F90 b/parteh/PRTAllometricCNPMod.F90 index eecace30b5..a8a0b826ff 100644 --- a/parteh/PRTAllometricCNPMod.F90 +++ b/parteh/PRTAllometricCNPMod.F90 @@ -258,6 +258,7 @@ module PRTAllometricCNPMod procedure :: DailyPRT => DailyPRTAllometricCNP procedure :: FastPRT => FastPRTAllometricCNP procedure :: GetNutrientTarget => GetNutrientTargetCNP + procedure :: CacheIndices => CacheIndicesAllometricCNP ! Extended functions specific to Allometric CNP procedure :: CNPPrioritizedReplacement @@ -367,9 +368,63 @@ subroutine InitPRTGlobalAllometricCNP() return end subroutine InitPRTGlobalAllometricCNP - ! ===================================================================================== + + subroutine CacheIndicesAllometricCNP(this) + + ! Cache variable indices for fast repeated access + ! Call this once after InitAllocate + + + class(cnp_allom_prt_vartypes) :: this + ! Carbon indices + this%i_leaf_c = prt_global%sp_organ_map(leaf_organ, carbon12_element) + this%i_fnrt_c = prt_global%sp_organ_map(fnrt_organ, carbon12_element) + this%i_sapw_c = prt_global%sp_organ_map(sapw_organ, carbon12_element) + this%i_store_c = prt_global%sp_organ_map(store_organ, carbon12_element) + this%i_struct_c = prt_global%sp_organ_map(struct_organ, carbon12_element) + this%i_repro_c = prt_global%sp_organ_map(repro_organ, carbon12_element) + + ! Carbon pointers + this%fnrt_c => this%variables(this%i_fnrt_c)%val(1) + this%sapw_c => this%variables(this%i_sapw_c)%val(1) + this%store_c => this%variables(this%i_store_c)%val(1) + this%struct_c => this%variables(this%i_struct_c)%val(1) + this%repro_c => this%variables(this%i_repro_c)%val(1) + this%leaf_c => this%variables(this%i_leaf_c)%val(:) + + ! Nitrogen indices + this%i_leaf_n = prt_global%sp_organ_map(leaf_organ, nitrogen_element) + this%i_fnrt_n = prt_global%sp_organ_map(fnrt_organ, nitrogen_element) + this%i_sapw_n = prt_global%sp_organ_map(sapw_organ, nitrogen_element) + this%i_store_n = prt_global%sp_organ_map(store_organ, nitrogen_element) + this%i_struct_n = prt_global%sp_organ_map(struct_organ, nitrogen_element) + + ! Nitrogen pointers + this%fnrt_n => this%variables(this%i_fnrt_n)%val(1) + this%sapw_n => this%variables(this%i_sapw_n)%val(1) + this%store_n => this%variables(this%i_store_n)%val(1) + this%struct_n => this%variables(this%i_struct_n)%val(1) + this%leaf_n => this%variables(this%i_leaf_n)%val(:) + + ! Phosphorus indices + this%i_leaf_p = prt_global%sp_organ_map(leaf_organ, phosphorus_element) + this%i_fnrt_p = prt_global%sp_organ_map(fnrt_organ, phosphorus_element) + this%i_sapw_p = prt_global%sp_organ_map(sapw_organ, phosphorus_element) + this%i_store_p = prt_global%sp_organ_map(store_organ, phosphorus_element) + this%i_struct_p = prt_global%sp_organ_map(struct_organ, phosphorus_element) + + ! Phosphorus pointers + this%fnrt_p => this%variables(this%i_fnrt_p)%val(1) + this%sapw_p => this%variables(this%i_sapw_p)%val(1) + this%store_p => this%variables(this%i_store_p)%val(1) + this%struct_p => this%variables(this%i_struct_p)%val(1) + this%leaf_p => this%variables(this%i_leaf_p)%val(:) + + + end subroutine CacheIndicesAllometricCNP + subroutine DailyPRTAllometricCNP(this,phase) diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index c5e46f783e..51163799b9 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -139,7 +139,8 @@ module PRTAllometricCarbonMod procedure :: DailyPRT => DailyPRTAllometricCarbon procedure :: FastPRT => FastPRTAllometricCarbon - + procedure :: CacheIndices => CacheIndicesAllometricCarbon + end type callom_prt_vartypes ! ------------------------------------------------------------------------------------ @@ -256,7 +257,32 @@ end subroutine InitPRTGlobalAllometricCarbon ! ===================================================================================== - + subroutine CacheIndicesAllometricCarbon(this) + + ! Cache variable indices for fast repeated access + ! Call this once after InitAllocate + + class(callom_prt_vartypes) :: this + + ! Carbon indices + this%i_leaf_c = prt_global%sp_organ_map(leaf_organ, carbon12_element) + this%i_fnrt_c = prt_global%sp_organ_map(fnrt_organ, carbon12_element) + this%i_sapw_c = prt_global%sp_organ_map(sapw_organ, carbon12_element) + this%i_store_c = prt_global%sp_organ_map(store_organ, carbon12_element) + this%i_struct_c = prt_global%sp_organ_map(struct_organ, carbon12_element) + this%i_repro_c = prt_global%sp_organ_map(repro_organ, carbon12_element) + + ! Carbon pointers + this%fnrt_c => this%variables(this%i_fnrt_c)%val(1) + this%sapw_c => this%variables(this%i_sapw_c)%val(1) + this%store_c => this%variables(this%i_store_c)%val(1) + this%struct_c => this%variables(this%i_struct_c)%val(1) + this%repro_c => this%variables(this%i_repro_c)%val(1) + this%leaf_c => this%variables(this%i_leaf_c)%val(:) + + + end subroutine CacheIndicesAllometricCarbon + subroutine DailyPRTAllometricCarbon(this,phase) ! ----------------------------------------------------------------------------------- diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 573a99cf3a..c1d3ac3015 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -178,7 +178,7 @@ module PRTGenericMod type, public :: prt_vartype - real(r8),pointer :: val(:) ! Instantaneous state variable [kg] + real(r8),pointer :: val(:) ! Instantaneous state variable [kg] real(r8),allocatable :: val0(:) ! State variable at the beginning ! of the control period [kg] real(r8),allocatable :: net_alloc(:) ! Net change due to allocation/transport [kg] @@ -230,6 +230,31 @@ module PRTGenericMod ! ------------------------------------------------------------------------------------- type, public :: prt_vartypes + + ! Performance optimization: Cached variable indices + integer :: i_leaf_c, i_fnrt_c, i_sapw_c, i_store_c, i_struct_c, i_repro_c + integer :: i_leaf_n, i_fnrt_n, i_sapw_n, i_store_n, i_struct_n + integer :: i_leaf_p, i_fnrt_p, i_sapw_p, i_store_p, i_struct_p + + real(r8), pointer :: fnrt_c => null() + real(r8), pointer :: sapw_c => null() + real(r8), pointer :: store_c => null() + real(r8), pointer :: struct_c => null() + real(r8), pointer :: repro_c => null() + + real(r8), pointer :: fnrt_n => null() + real(r8), pointer :: sapw_n => null() + real(r8), pointer :: store_n => null() + real(r8), pointer :: struct_n => null() + + real(r8), pointer :: fnrt_p => null() + real(r8), pointer :: sapw_p => null() + real(r8), pointer :: store_p => null() + real(r8), pointer :: struct_p => null() + + real(r8), pointer :: leaf_c(:) => null() + real(r8), pointer :: leaf_n(:) => null() + real(r8), pointer :: leaf_p(:) => null() type(prt_vartype),allocatable :: variables(:) ! The state variables and fluxes type(prt_bctype), allocatable :: bc_inout(:) ! These boundaries may be changed @@ -246,7 +271,8 @@ module PRTGenericMod procedure :: FastPRT => FastPRTBase procedure :: DamageRecovery => DamageRecoveryBase procedure :: GetNutrientTarget => GetNutrientTargetBase - + procedure :: CacheIndices => CacheIndicesBase + ! These are generic functions that should work on all hypotheses procedure, non_overridable :: InitAllocate @@ -269,6 +295,15 @@ module PRTGenericMod procedure, non_overridable :: WeightedFusePRTVartypes procedure, non_overridable :: CopyPRTVartypes + procedure, non_overridable :: GetLeafC + procedure, non_overridable :: GetFnrtC + procedure, non_overridable :: GetSapwC + procedure, non_overridable :: GetStoreC + procedure, non_overridable :: GetStructC + procedure, non_overridable :: GetReproC + procedure, non_overridable :: GetLeafN + procedure, non_overridable :: GetSapwN + procedure, non_overridable :: GetFnrtN procedure :: AgeLeaves ! This routine may be used generically ! but also leaving the door open for over-rides @@ -554,6 +589,8 @@ subroutine InitAllocate(this) end do + ! Cache variable indices for fast access (performance optimization) + call this%CacheIndices() return end subroutine InitAllocate @@ -1057,6 +1094,69 @@ function GetState(this, organ_id, element_id, position_id) result(state_val) return end function GetState + ! =================================================================================== + + + + ! Fast carbon accessors (no map lookup, direct index access) + + function GetLeafC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = sum(this%variables(this%i_leaf_c)%val(:)) + end function GetLeafC + + function GetFnrtC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_fnrt_c)%val(1) + end function GetFnrtC + + function GetSapwC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_sapw_c)%val(1) + end function GetSapwC + + function GetStoreC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_store_c)%val(1) + end function GetStoreC + + function GetStructC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_struct_c)%val(1) + end function GetStructC + + function GetReproC(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_repro_c)%val(1) + end function GetReproC + + ! Fast nitrogen accessors + + function GetLeafN(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_leaf_n)%val(1) + end function GetLeafN + + function GetSapwN(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_sapw_n)%val(1) + end function GetSapwN + + function GetFnrtN(this) result(val) + class(prt_vartypes), intent(in) :: this + real(r8) :: val + val = this%variables(this%i_fnrt_n)%val(1) + end function GetFnrtN + + ! ==================================================================================== subroutine GetBiomass(this, element_id, & @@ -1252,6 +1352,14 @@ function GetCoordVal(this, organ_id, element_id ) result(prt_val) call endrun(msg=errMsg(sourcefile, __LINE__)) end function GetCoordVal + + subroutine CacheIndicesBase(this) + + class(prt_vartypes) :: this + write(fates_log(),*)'CacheIndices must be extended' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + end subroutine CacheIndicesBase ! ==================================================================================== diff --git a/radiation/FatesNormanRadMod.F90 b/radiation/FatesNormanRadMod.F90 index c16635180a..0085f2e706 100644 --- a/radiation/FatesNormanRadMod.F90 +++ b/radiation/FatesNormanRadMod.F90 @@ -191,7 +191,7 @@ subroutine PatchNormanRadiation (currentPatch, & tau_layer(:,:,:,:)=0.0_r8 f_abs(:,:,:,:)=0.0_r8 f_abs_leaf(:,:,:,:)=0._r8 - do L = 1,currentPatch%ncl_p + do L = 1,currentPatch%ncl do ft = 1,numpft currentPatch%canopy_mask(L,ft) = 0 do iv = 1, currentPatch%nrad(L,ft) @@ -257,7 +257,7 @@ subroutine PatchNormanRadiation (currentPatch, & ! ------------------------------------------------------------------------------ ftweight(:,:,:) = 0._r8 - do L = 1,currentPatch%NCL_p + do L = 1,currentPatch%ncl do ft = 1,numpft do iv = 1, currentPatch%nrad(L,ft) !this is already corrected for area in CLAP @@ -275,7 +275,7 @@ subroutine PatchNormanRadiation (currentPatch, & endif end if - do L = 1,currentPatch%NCL_p !start at the top canopy layer (1 is the top layer.) + do L = 1,currentPatch%ncl !start at the top canopy layer (1 is the top layer.) weighted_dir_tr(L) = 0.0_r8 weighted_fsun(L) = 0._r8 @@ -422,7 +422,7 @@ subroutine PatchNormanRadiation (currentPatch, & end do !L - do L = currentPatch%NCL_p,1, -1 !start at the bottom and work up. + do L = currentPatch%ncl,1, -1 !start at the bottom and work up. do ft = 1,numpft if (currentPatch%canopy_mask(L,ft) == 1)then @@ -449,7 +449,7 @@ subroutine PatchNormanRadiation (currentPatch, & !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! ! Soil diffuse reflectance (ratio of down to up radiation). iv = currentPatch%nrad(L,ft) + 1 - if (L == currentPatch%NCL_p)then !nearest the soil + if (L == currentPatch%ncl)then !nearest the soil dif_ratio(L,ft,iv,ib) = currentPatch%gnd_alb_dif(ib) !bc_in(s)%albgr_dif_rb(ib) else dif_ratio(L,ft,iv,ib) = weighted_dif_ratio(L+1,ib) @@ -481,7 +481,7 @@ subroutine PatchNormanRadiation (currentPatch, & Dif_dn(:,:,:) = 0.00_r8 Dif_up(:,:,:) = 0.00_r8 - do L = 1, currentPatch%NCL_p !work down from the top of the canopy. + do L = 1, currentPatch%ncl !work down from the top of the canopy. weighted_dif_down(L) = 0._r8 do ft = 1, numpft if (currentPatch%canopy_mask(L,ft) == 1)then @@ -530,20 +530,20 @@ subroutine PatchNormanRadiation (currentPatch, & !instance where the first layer ftweight is used a proxy for the whole column. FTWA endif !present end do !ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is the the (incomplete) understorey? + if (L == currentPatch%ncl.and.currentPatch%ncl > 1)then !is the the (incomplete) understorey? !Add on the radiation going through the canopy gaps. weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1)*(1.0-sum(ftweight(L,:,1))) !instance where the first layer ftweight is used a proxy for the whole column. FTWA endif end do !L - do L = currentPatch%NCL_p,1 ,-1 !work up from the bottom. + do L = currentPatch%ncl,1 ,-1 !work up from the bottom. weighted_dif_up(L) = 0._r8 do ft = 1, numpft if (currentPatch%canopy_mask(L,ft) == 1)then !Bounce diffuse radiation off soil surface. iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then !is this the bottom layer ? + if (L==currentPatch%ncl)then !is this the bottom layer ? Dif_up(L,ft,iv) = currentPatch%gnd_alb_dif(ib) * Dif_dn(L,ft,iv) else Dif_up(L,ft,iv) = weighted_dif_up(L+1) @@ -570,7 +570,7 @@ subroutine PatchNormanRadiation (currentPatch, & !instance where the first layer ftweight is used a proxy for the whole column. FTWA endif !present end do !ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + if (L == currentPatch%ncl.and.currentPatch%ncl > 1)then !is this the (incomplete) understorey? !Add on the radiation coming up through the canopy gaps. !diffuse to diffuse weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & @@ -594,7 +594,7 @@ subroutine PatchNormanRadiation (currentPatch, & iter = iter + 1 irep = 0 - do L = 1,currentPatch%NCL_p !working from the top down + do L = 1,currentPatch%ncl !working from the top down weighted_dif_down(L) = 0._r8 do ft =1,numpft if (currentPatch%canopy_mask(L,ft) == 1)then @@ -651,19 +651,19 @@ subroutine PatchNormanRadiation (currentPatch, & endif !present end do!ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + if (L == currentPatch%ncl.and.currentPatch%ncl > 1)then !is this the (incomplete) understorey? weighted_dif_down(L) = weighted_dif_down(L) + weighted_dif_down(L-1) * & (1.0_r8-sum(ftweight(L,1:numpft,1))) end if end do ! do L loop - do L = 1, currentPatch%NCL_p ! working from the top down. + do L = 1, currentPatch%ncl ! working from the top down. weighted_dif_up(L) = 0._r8 do ft =1,numpft if (currentPatch%canopy_mask(L,ft) == 1)then ! Upward diffuse flux at soil or from lower canopy (forward diffuse and unscattered direct beam) iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then !In the bottom canopy layer, reflect off the soil + if (L==currentPatch%ncl)then !In the bottom canopy layer, reflect off the soil Dif_up(L,ft,iv) = Dif_dn(L,ft,iv) * currentPatch%gnd_alb_dif(ib) + & forc_dir(radtype) * tr_dir_z(L,ft,iv) * currentPatch%gnd_alb_dir(ib) else !In the other canopy layers, reflect off the underlying vegetation. @@ -695,7 +695,7 @@ subroutine PatchNormanRadiation (currentPatch, & end if !present end do!ft - if (L == currentPatch%NCL_p.and.currentPatch%NCL_p > 1)then !is this the (incomplete) understorey? + if (L == currentPatch%ncl.and.currentPatch%ncl > 1)then !is this the (incomplete) understorey? !Add on the radiation coming up through the canopy gaps. weighted_dif_up(L) = weighted_dif_up(L) +(1.0_r8-sum(ftweight(L,1:numpft,1))) * & weighted_dif_down(L-1) * currentPatch%gnd_alb_dif(ib) @@ -709,7 +709,7 @@ subroutine PatchNormanRadiation (currentPatch, & tr_soili = 0._r8 tr_soild = 0._r8 - do L = 1, currentPatch%NCL_p !working from the top down. + do L = 1, currentPatch%ncl !working from the top down. abs_dir_z(:,:) = 0._r8 abs_dif_z(:,:) = 0._r8 do ft =1,numpft @@ -728,7 +728,7 @@ subroutine PatchNormanRadiation (currentPatch, & end do ! Absorbed direct beam and diffuse do soil - if (L == currentPatch%NCL_p)then + if (L == currentPatch%ncl)then iv = currentPatch%nrad(L,ft) + 1 Abs_dif_z(ft,iv) = ftweight(L,ft,1)*Dif_dn(L,ft,iv) * (1.0_r8 - currentPatch%gnd_alb_dif(ib) ) Abs_dir_z(ft,iv) = ftweight(L,ft,1)*forc_dir(radtype) * & @@ -772,7 +772,7 @@ subroutine PatchNormanRadiation (currentPatch, & !==============================================================================! ! Solar radiation absorbed by ground iv = currentPatch%nrad(L,ft) + 1 - if (L==currentPatch%NCL_p)then + if (L==currentPatch%ncl)then abs_rad(ib) = abs_rad(ib) + (Abs_dir_z(ft,iv) + Abs_dif_z(ft,iv)) end if ! Solar radiation absorbed by vegetation and sunlit/shaded leaves @@ -819,7 +819,7 @@ subroutine PatchNormanRadiation (currentPatch, & !radiation absorbed from fluxes through unfilled part of lower canopy. - if (currentPatch%NCL_p > 1.and.L == currentPatch%NCL_p)then + if (currentPatch%ncl > 1.and.L == currentPatch%ncl)then abs_rad(ib) = abs_rad(ib) + weighted_dif_down(L-1) * & (1.0_r8-sum(ftweight(L,1:numpft,1)))*(1.0_r8-currentPatch%gnd_alb_dif(ib) ) abs_rad(ib) = abs_rad(ib) + forc_dir(radtype) * weighted_dir_tr(L-1) * & @@ -857,7 +857,7 @@ subroutine PatchNormanRadiation (currentPatch, & if ( abs(error) > 0.0001)then write(fates_log(),*)'dir ground absorption error',error,currentPatch%sabs_dir(ib), & currentPatch%tr_soil_dir(ib)* & - (1.0_r8-currentPatch%gnd_alb_dir(ib) ),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) + (1.0_r8-currentPatch%gnd_alb_dir(ib) ),currentPatch%ncl,ib,sum(ftweight(1,1:numpft,1)) write(fates_log(),*) 'albedos',currentPatch%sabs_dir(ib) ,currentPatch%tr_soil_dir(ib), & (1.0_r8-currentPatch%gnd_alb_dir(ib) ) do ft =1,numpft @@ -873,7 +873,7 @@ subroutine PatchNormanRadiation (currentPatch, & (1.0_r8-currentPatch%gnd_alb_dif(ib) ))) > 0.0001_r8)then write(fates_log(),*)'dif ground absorption error',currentPatch%sabs_dif(ib) , & (currentPatch%tr_soil_dif(ib)* & - (1.0_r8-currentPatch%gnd_alb_dif(ib) )),currentPatch%NCL_p,ib,sum(ftweight(1,1:numpft,1)) + (1.0_r8-currentPatch%gnd_alb_dif(ib) )),currentPatch%ncl,ib,sum(ftweight(1,1:numpft,1)) endif end if endif @@ -896,7 +896,7 @@ subroutine PatchNormanRadiation (currentPatch, & endif lai_reduction(:) = 0.0_r8 - do L = 1, currentPatch%NCL_p + do L = 1, currentPatch%ncl do ft =1,numpft if (currentPatch%canopy_mask(L,ft) == 1)then do iv = 1, currentPatch%nrad(L,ft) @@ -926,8 +926,8 @@ subroutine PatchNormanRadiation (currentPatch, & write(fates_log(),*) 'Large Dir Radn consvn error',error ,ib write(fates_log(),*) 'diags', albd_parb_out(ib), ftdd_parb_out(ib), & ftid_parb_out(ib), fabd_parb_out(ib) - write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl,1:numpft,1:diag_nlevleaf) + write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl,1:numpft,1:diag_nlevleaf) write(fates_log(),*) 'ftweight',ftweight(1,1:numpft,1:diag_nlevleaf) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) @@ -945,10 +945,10 @@ subroutine PatchNormanRadiation (currentPatch, & write(fates_log(),*) 'lg Dif Radn consvn error',error ,ib write(fates_log(),*) 'diags', albi_parb_out(ib), ftii_parb_out(ib), & fabi_parb_out(ib) - !write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - !write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - !write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) - !write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl_p,1:numpft,1:diag_nlevleaf) + !write(fates_log(),*) 'lai_change',lai_change(currentpatch%ncl,1:numpft,1:diag_nlevleaf) + !write(fates_log(),*) 'elai',currentpatch%elai_profile(currentpatch%ncl,1:numpft,1:diag_nlevleaf) + !write(fates_log(),*) 'esai',currentpatch%esai_profile(currentpatch%ncl,1:numpft,1:diag_nlevleaf) + !write(fates_log(),*) 'ftweight',ftweight(currentpatch%ncl,1:numpft,1:diag_nlevleaf) write(fates_log(),*) 'cp',currentPatch%area, currentPatch%patchno write(fates_log(),*) 'ground albedo diffuse (ib)', currentPatch%gnd_alb_dir(ib) !write(fates_log(),*) 'rhol',rhol(1:numpft,:) diff --git a/radiation/FatesRadiationDriveMod.F90 b/radiation/FatesRadiationDriveMod.F90 index a401388f7e..f63adc5fa9 100644 --- a/radiation/FatesRadiationDriveMod.F90 +++ b/radiation/FatesRadiationDriveMod.F90 @@ -42,20 +42,21 @@ module FatesRadiationDriveMod use FatesGlobals, only : endrun => fates_endrun use EDPftvarcon, only : EDPftvarcon_inst use FatesNormanRadMod, only : PatchNormanRadiation - + ! CIME globals use shr_log_mod , only : errMsg => shr_log_errMsg implicit none private + public :: FatesNormalizedPatchRadiation public :: FatesNormalizedCanopyRadiation ! Surface albedo and two-stream fluxes public :: FatesSunShadeFracs - logical :: debug = .false. ! for debugging this module + logical,parameter :: debug = .false. ! for debugging this module character(len=*), parameter, private :: sourcefile = & __FILE__ - + contains subroutine FatesNormalizedCanopyRadiation(sites, bc_in, bc_out ) @@ -82,14 +83,7 @@ subroutine FatesNormalizedCanopyRadiation(sites, bc_in, bc_out ) integer :: s ! site loop counter integer :: nsites ! number of sites integer :: ifp ! patch loop counter - integer :: ib ! radiation broad band counter - integer :: cl, iv, icol, ft ! indices for canopy layer,leaf layer, - ! rad column and functional type - integer :: nv ! number of veg layers - real(r8) :: area_frac ! area fraction for layer of interest - real(r8) :: vai_top ! integrated (top-down) vegetation area - ! index at lop of layer - real(r8) :: vai ! total VAI of the scattering element + type(fates_patch_type), pointer :: currentPatch ! patch pointer !----------------------------------------------------------------------- @@ -99,7 +93,7 @@ subroutine FatesNormalizedCanopyRadiation(sites, bc_in, bc_out ) ! re-examined int he future ! RGK,2016-08-06: FATES is still incompatible with VOC emission module ! ------------------------------------------------------------------------------- - + nsites = size(sites,dim=1) do s = 1, nsites @@ -111,127 +105,150 @@ subroutine FatesNormalizedCanopyRadiation(sites, bc_in, bc_out ) do while (associated(currentpatch)) ifp = currentpatch%patchno - + if_bareground: if(currentpatch%nocomp_pft_label.ne.nocomp_bareground)then - - ! Initialize output boundary conditions with trivial assumption - ! This matches CLM/ELM - ! Albedo is perfect reflector, no flux into or through canopy - bc_out(s)%albd_parb(ifp,:) = 1._r8 - bc_out(s)%albi_parb(ifp,:) = 1._r8 - bc_out(s)%fabi_parb(ifp,:) = 0._r8 - bc_out(s)%fabd_parb(ifp,:) = 0._r8 - bc_out(s)%ftdd_parb(ifp,:) = 0._r8 - bc_out(s)%ftid_parb(ifp,:) = 0._r8 - bc_out(s)%ftii_parb(ifp,:) = 0._r8 - - ! Zero diagnostics - currentPatch%f_sun (:,:,:) = 0._r8 - currentPatch%fabd_sun_z (:,:,:) = 0._r8 - currentPatch%fabd_sha_z (:,:,:) = 0._r8 - currentPatch%fabi_sun_z (:,:,:) = 0._r8 - currentPatch%fabi_sha_z (:,:,:) = 0._r8 - currentPatch%fabd (:) = 0._r8 - currentPatch%fabi (:) = 0._r8 - currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:) = 0._r8 - currentPatch%nrmlzd_parprof_pft_dif_z(:,:,:) = 0._r8 - currentPatch%gnd_alb_dif(1:num_swb) = bc_in(s)%albgr_dif_rb(1:num_swb) - currentPatch%gnd_alb_dir(1:num_swb) = bc_in(s)%albgr_dir_rb(1:num_swb) - currentPatch%fcansno = bc_in(s)%fcansno_pa(ifp) - currentPatch%rad_error(:) = hlm_hio_ignore_val - - if_zenith_flag: if( bc_in(s)%coszen>0._r8 )then - - select case(hlm_radiation_model) - case(norman_solver) - - call PatchNormanRadiation (currentPatch, & - bc_in(s)%coszen, & - bc_out(s)%albd_parb(ifp,:), & ! Surface Albedo direct - bc_out(s)%albi_parb(ifp,:), & ! Surface Albedo (indirect) diffuse - bc_out(s)%fabd_parb(ifp,:), & ! Fraction direct absorbed by canopy per unit incident - bc_out(s)%fabi_parb(ifp,:), & ! Fraction diffuse absorbed by canopy per unit incident - bc_out(s)%ftdd_parb(ifp,:), & ! Down direct flux below canopy per unit direct at top - bc_out(s)%ftid_parb(ifp,:), & ! Down diffuse flux below canopy per unit direct at top - bc_out(s)%ftii_parb(ifp,:)) ! Down diffuse flux below canopy per unit diffuse at top - - case(twostr_solver) - - associate( twostr => currentPatch%twostr) - - call twostr%CanopyPrep(currentPatch%fcansno) - call twostr%ZenithPrep(sites(s)%coszen) - - do ib = 1,num_swb - - twostr%band(ib)%albedo_grnd_diff = currentPatch%gnd_alb_dif(ib) - twostr%band(ib)%albedo_grnd_beam = currentPatch%gnd_alb_dir(ib) - - call twostr%Solve(ib, & ! in - normalized_upper_boundary, & ! in - 1.0_r8,1.0_r8, & ! in - sites(s)%taulambda_2str, & ! inout (scratch) - sites(s)%omega_2str, & ! inout (scratch) - sites(s)%ipiv_2str, & ! inout (scratch) - bc_out(s)%albd_parb(ifp,ib), & ! out - bc_out(s)%albi_parb(ifp,ib), & ! out - currentPatch%rad_error(ib), & ! out - bc_out(s)%fabd_parb(ifp,ib), & ! out - bc_out(s)%fabi_parb(ifp,ib), & ! out - bc_out(s)%ftdd_parb(ifp,ib), & ! out - bc_out(s)%ftid_parb(ifp,ib), & ! out - bc_out(s)%ftii_parb(ifp,ib)) - - if(debug) then - currentPatch%twostr%band(ib)%Rbeam_atm = 1._r8 - currentPatch%twostr%band(ib)%Rdiff_atm = 1._r8 - call CheckPatchRadiationBalance(currentPatch, sites(s)%snow_depth, & - ib, bc_out(s)%fabd_parb(ifp,ib),bc_out(s)%fabi_parb(ifp,ib)) - currentPatch%twostr%band(ib)%Rbeam_atm = fates_unset_r8 - currentPatch%twostr%band(ib)%Rdiff_atm = fates_unset_r8 - - if(bc_out(s)%fabi_parb(ifp,ib)>1.0 .or. bc_out(s)%fabd_parb(ifp,ib)>1.0)then - write(fates_log(),*) 'absorbed fraction > 1.0?' - write(fates_log(),*) ifp,ib,bc_out(s)%fabi_parb(ifp,ib),bc_out(s)%fabd_parb(ifp,ib) - call twostr%Dump(ib,lat=sites(s)%lat,lon=sites(s)%lon) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end if - end do - ! Fill in the diagnostic arrays for normalized radiation profiles - do_cl: do cl = 1,twostr%n_lyr - do_icol: do icol = 1,twostr%n_col(cl) - ft = twostr%scelg(cl,icol)%pft - if_notair: if (ft>0) then - area_frac = twostr%scelg(cl,icol)%area - vai = twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai - nv = GetNVegLayers(vai) - do iv = 1, nv - vai_top = dlower_vai(iv) - currentPatch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) = currentPatch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) + & - area_frac*twostr%GetRb(cl,icol,ivis,vai_top) - currentPatch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) = currentPatch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) + & - area_frac*twostr%GetRdDn(cl,icol,ivis,vai_top) + & - area_frac*twostr%GetRdUp(cl,icol,ivis,vai_top) - end do - end if if_notair - end do do_icol - end do do_cl - - end associate - end select - endif if_zenith_flag + call FatesNormalizedPatchRadiation(ifp,sites(s),bc_in(s),bc_out(s)) + end if if_bareground currentPatch => currentPatch%younger end do end do - + return end subroutine FatesNormalizedCanopyRadiation ! ====================================================================================== + subroutine FatesNormalizedPatchRadiation(ifp,site,bc_in,bc_out) + + type(ed_site_type), intent(inout) :: site + type(bc_in_type), intent(in) :: bc_in + type(bc_out_type), intent(inout) :: bc_out + integer,intent(in) :: ifp ! patch loop counter + + type(fates_patch_type), pointer :: patch + integer :: ib ! radiation broad band counter + integer :: cl, iv, icol, ft ! indices for canopy layer,leaf layer, + ! rad column and functional type + integer :: nv ! number of veg layers + real(r8) :: area_frac ! area fraction for layer of interest + real(r8) :: vai_top ! integrated (top-down) vegetation area + ! index at lop of layer + real(r8) :: vai ! total VAI of the scattering element + + patch => site%pa_vec(ifp)%p + + ! Initialize output boundary conditions with trivial assumption + ! This matches CLM/ELM + ! Albedo is perfect reflector, no flux into or through canopy + bc_out%albd_parb(ifp,:) = 1._r8 + bc_out%albi_parb(ifp,:) = 1._r8 + bc_out%fabi_parb(ifp,:) = 0._r8 + bc_out%fabd_parb(ifp,:) = 0._r8 + bc_out%ftdd_parb(ifp,:) = 0._r8 + bc_out%ftid_parb(ifp,:) = 0._r8 + bc_out%ftii_parb(ifp,:) = 0._r8 + + ! Zero diagnostics + patch%f_sun (:,:,:) = 0._r8 + patch%fabd_sun_z (:,:,:) = 0._r8 + patch%fabd_sha_z (:,:,:) = 0._r8 + patch%fabi_sun_z (:,:,:) = 0._r8 + patch%fabi_sha_z (:,:,:) = 0._r8 + patch%fabd (:) = 0._r8 + patch%fabi (:) = 0._r8 + patch%nrmlzd_parprof_pft_dir_z(:,:,:) = 0._r8 + patch%nrmlzd_parprof_pft_dif_z(:,:,:) = 0._r8 + patch%gnd_alb_dif(1:num_swb) = bc_in%albgr_dif_rb(1:num_swb) + patch%gnd_alb_dir(1:num_swb) = bc_in%albgr_dir_rb(1:num_swb) + patch%fcansno = bc_in%fcansno_pa(ifp) + patch%rad_error(:) = hlm_hio_ignore_val + + if_zenith_flag: if( bc_in%coszen>0._r8 )then + + select case(hlm_radiation_model) + case(norman_solver) + + call PatchNormanRadiation (patch, & + bc_in%coszen, & + bc_out%albd_parb(ifp,:), & ! Surface Albedo direct + bc_out%albi_parb(ifp,:), & ! Surface Albedo (indirect) diffuse + bc_out%fabd_parb(ifp,:), & ! Fraction direct absorbed by canopy per unit incident + bc_out%fabi_parb(ifp,:), & ! Fraction diffuse absorbed by canopy per unit incident + bc_out%ftdd_parb(ifp,:), & ! Down direct flux below canopy per unit direct at top + bc_out%ftid_parb(ifp,:), & ! Down diffuse flux below canopy per unit direct at top + bc_out%ftii_parb(ifp,:)) ! Down diffuse flux below canopy per unit diffuse at top + + case(twostr_solver) + + associate( twostr => patch%twostr) + + call twostr%CanopyPrep(patch%fcansno) + call twostr%ZenithPrep(site%coszen) + + do ib = 1,num_swb + + twostr%band(ib)%albedo_grnd_diff = patch%gnd_alb_dif(ib) + twostr%band(ib)%albedo_grnd_beam = patch%gnd_alb_dir(ib) + + call twostr%Solve(ib, & ! in + normalized_upper_boundary, & ! in + 1.0_r8,1.0_r8, & ! in + bc_out%albd_parb(ifp,ib), & ! out + bc_out%albi_parb(ifp,ib), & ! out + patch%rad_error(ib), & ! out + bc_out%fabd_parb(ifp,ib), & ! out + bc_out%fabi_parb(ifp,ib), & ! out + bc_out%ftdd_parb(ifp,ib), & ! out + bc_out%ftid_parb(ifp,ib), & ! out + bc_out%ftii_parb(ifp,ib)) + + if(debug) then + patch%twostr%band(ib)%Rbeam_atm = 1._r8 + patch%twostr%band(ib)%Rdiff_atm = 1._r8 + call CheckPatchRadiationBalance(patch, site%snow_depth, & + ib, bc_out%fabd_parb(ifp,ib),bc_out%fabi_parb(ifp,ib)) + patch%twostr%band(ib)%Rbeam_atm = fates_unset_r8 + patch%twostr%band(ib)%Rdiff_atm = fates_unset_r8 + + if(bc_out%fabi_parb(ifp,ib)>1.0 .or. bc_out%fabd_parb(ifp,ib)>1.0)then + write(fates_log(),*) 'absorbed fraction > 1.0?' + write(fates_log(),*) ifp,ib,bc_out%fabi_parb(ifp,ib),bc_out%fabd_parb(ifp,ib) + call twostr%Dump(ib,lat=site%lat,lon=site%lon) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end if + end do + + ! Fill in the diagnostic arrays for normalized radiation profiles + do_cl: do cl = 1,twostr%n_lyr + do_icol: do icol = 1,twostr%n_col(cl) + ft = twostr%scelg(cl,icol)%pft + if_notair: if (ft>0) then + area_frac = twostr%scelg(cl,icol)%area + vai = twostr%scelg(cl,icol)%sai+twostr%scelg(cl,icol)%lai + nv = GetNVegLayers(vai) + do iv = 1, nv + vai_top = dlower_vai(iv) + patch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) = patch%nrmlzd_parprof_pft_dir_z(cl,ft,iv) + & + area_frac*twostr%GetRb(cl,icol,ivis,vai_top) + patch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) = patch%nrmlzd_parprof_pft_dif_z(cl,ft,iv) + & + area_frac*twostr%GetRdDn(cl,icol,ivis,vai_top) + & + area_frac*twostr%GetRdUp(cl,icol,ivis,vai_top) + end do + end if if_notair + end do do_icol + end do do_cl + + end associate + end select + endif if_zenith_flag + + end subroutine FatesNormalizedPatchRadiation + + ! ========================================================================================= + subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) implicit none @@ -262,18 +279,18 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) real(r8) :: vai logical :: call_fail type(fates_patch_type), pointer :: fpatch ! patch pointer for failure reporting - + do s = 1,nsites cpatch => sites(s)%oldest_patch do while (associated(cpatch)) ifp = cpatch%patchno - + if_bareground:if(cpatch%nocomp_pft_label.ne.nocomp_bareground)then !only for veg patches ! do not do albedo calculations for bare ground patch in SP mode - + ! Initialize diagnostics cpatch%ed_parsun_z(:,:,:) = 0._r8 cpatch%ed_parsha_z(:,:,:) = 0._r8 @@ -284,14 +301,14 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) sunlai = 0._r8 shalai = 0._r8 - + ! Loop over patches to calculate laisun_z and laisha_z for each layer. ! Derive canopy laisun, laisha, and fsun from layer sums. ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value. ! cpatch%f_sun is calculated in the surface_albedo routine... - - do cl = 1, cpatch%ncl_p + + do cl = 1, cpatch%ncl do ft = 1,numpft do iv = 1,cpatch%nrad(cl,ft) cpatch%ed_laisun_z(cl,ft,iv) = cpatch%elai_profile(cl,ft,iv) * & @@ -310,41 +327,41 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) else bc_out(s)%fsun_pa(ifp) = 0._r8 endif - + if(debug)then if(bc_out(s)%fsun_pa(ifp) > 1._r8)then write(fates_log(),*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & sunlai,shalai endif end if - + elai = calc_areaindex(cpatch,'elai') - + bc_out(s)%laisun_pa(ifp) = elai*bc_out(s)%fsun_pa(ifp) bc_out(s)%laisha_pa(ifp) = elai*(1.0_r8-bc_out(s)%fsun_pa(ifp)) - + ! Absorbed PAR profile through canopy ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo ! are canopy integrated so that layer values equal big leaf values. - - do cl = 1, cpatch%ncl_p + + do cl = 1, cpatch%ncl do ft = 1,numpft do iv = 1, cpatch%nrad(cl,ft) - + cpatch%ed_parsun_z(cl,ft,iv) = & bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sun_z(cl,ft,iv) + & bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sun_z(cl,ft,iv) - + cpatch%ed_parsha_z(cl,ft,iv) = & bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sha_z(cl,ft,iv) + & bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sha_z(cl,ft,iv) - + end do !iv end do !ft end do !cl - - - + + + else ! if_norm_twostr ! If there is no sun out, we have a trivial solution @@ -354,7 +371,7 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) bc_out(s)%laisun_pa(ifp) = 0._r8 bc_out(s)%laisha_pa(ifp) = calc_areaindex(cpatch,'elai') bc_out(s)%fsun_pa(ifp) = 0._r8 - + else ! Two-stream @@ -363,20 +380,20 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) cpatch%twostr%band(ib)%Rbeam_atm = bc_in(s)%solad_parb(ifp,ib) cpatch%twostr%band(ib)%Rdiff_atm = bc_in(s)%solai_parb(ifp,ib) end do - + area_vlpfcl(:,:,:) = 0._r8 cpatch%f_sun(:,:,:) = 0._r8 - + call FatesPatchFSun(sites(s),cpatch, & bc_out(s)%fsun_pa(ifp), & bc_out(s)%laisun_pa(ifp), & bc_out(s)%laisha_pa(ifp)) - + associate(twostr => cpatch%twostr) - + do_cl: do cl = 1,twostr%n_lyr do_icol: do icol = 1,twostr%n_col(cl) - + ft = twostr%scelg(cl,icol)%pft if_notair: if (ft>0) then area_frac = twostr%scelg(cl,icol)%area @@ -385,7 +402,7 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) nv = GetNVegLayers(vai) do iv = 1, nv - + vai_top = dlower_vai(iv) if(iv == nv) then @@ -393,7 +410,7 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) else vai_bot = dlower_vai(iv+1) end if - + call twostr%GetAbsRad(cl,icol,ipar,vai_top,vai_bot, & Rb_abs,Rd_abs,Rd_abs_leaf,Rb_abs_leaf,R_abs_stem,R_abs_snow,leaf_sun_frac,call_fail) @@ -407,19 +424,19 @@ subroutine FatesSunShadeFracs(nsites, sites,bc_in,bc_out) call twostr%Dump(ipar,lat=sites(s)%lat,lon=sites(s)%lon) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + cpatch%f_sun(cl,ft,iv) = cpatch%f_sun(cl,ft,iv) + & area_frac*leaf_sun_frac cpatch%ed_parsun_z(cl,ft,iv) = cpatch%ed_parsun_z(cl,ft,iv) + & area_frac*(rd_abs_leaf*leaf_sun_frac + rb_abs_leaf) cpatch%ed_parsha_z(cl,ft,iv) = cpatch%ed_parsha_z(cl,ft,iv) + & area_frac*rd_abs_leaf*(1._r8-leaf_sun_frac) - + area_vlpfcl(iv,ft,cl) = area_vlpfcl(iv,ft,cl) + area_frac end do end if if_notair end do do_icol - + do ft = 1,numpft do_iv: do iv = 1,cpatch%nleaf(cl,ft) if(area_vlpfcl(iv,ft,cl) cpatch%younger enddo diff --git a/radiation/FatesTwoStreamUtilsMod.F90 b/radiation/FatesTwoStreamUtilsMod.F90 index f8d496c8b1..a0c93bdd7b 100644 --- a/radiation/FatesTwoStreamUtilsMod.F90 +++ b/radiation/FatesTwoStreamUtilsMod.F90 @@ -77,8 +77,6 @@ subroutine FatesConstructRadElements(site) ! and elements as well (to preserve mass and volume). integer :: max_elements ! Maximum number of scattering elements on the site - integer :: n_scr ! The size of the scratch arrays - logical :: allocate_scratch ! Whether to re-allocate the scratch arrays integer :: icolmax ! Column index for each layer with largest area footprint real(r8) :: areamax ! The area footprint of the largest column @@ -162,7 +160,7 @@ subroutine FatesConstructRadElements(site) end if ! Add the air element if the canopy area is not filled - do ican = 1,patch%ncl_p + do ican = 1,patch%ncl if( (1._r8-canopy_frac(ican))>area_err_thresh ) then n_col(ican) = n_col(ican) + 1 end if @@ -176,25 +174,25 @@ subroutine FatesConstructRadElements(site) ! ------------------------------------------------------------------------------------------- maxcol = 0 - do ican = 1,patch%ncl_p + do ican = 1,patch%ncl if (n_col(ican)>maxcol) maxcol=n_col(ican) end do if(.not.associated(twostr%scelg)) then - call twostr%AllocInitTwoStream((/ivis,inir/),patch%ncl_p,maxcol+2) + call twostr%AllocInitTwoStream((/ivis,inir/),patch%ncl,maxcol+2) else if(ubound(twostr%scelg,2) < maxcol .or. & ubound(twostr%scelg,2) > (maxcol+4) .or. & - ubound(twostr%scelg,1) < patch%ncl_p ) then + ubound(twostr%scelg,1) < patch%ncl ) then call twostr%DeallocTwoStream() ! Add a little more space than necessary so ! we don't have to keep allocating/deallocating - call twostr%AllocInitTwoStream((/ivis,inir/),patch%ncl_p,maxcol+2) + call twostr%AllocInitTwoStream((/ivis,inir/),patch%ncl,maxcol+2) end if @@ -266,7 +264,7 @@ subroutine FatesConstructRadElements(site) enddo - do ican = 1,patch%ncl_p + do ican = 1,patch%ncl ! If the canopy is not full, add an air element if( (1._r8-canopy_frac(ican))>area_err_thresh ) then @@ -332,12 +330,12 @@ subroutine FatesConstructRadElements(site) end do - twostr%n_col(1:patch%ncl_p) = n_col(1:patch%ncl_p) + twostr%n_col(1:patch%ncl) = n_col(1:patch%ncl) ! Set up some non-element parameters ! ------------------------------------------------------------------------------------------- - twostr%n_lyr = patch%ncl_p ! Number of layers + twostr%n_lyr = patch%ncl ! Number of layers call twostr%GetNSCel() ! Total number of elements @@ -352,33 +350,6 @@ subroutine FatesConstructRadElements(site) patch => patch%younger end do - ! Re-evaluate the scratch space used for solving two-stream radiation - ! The scratch space needs to be 2x the number of computational elements - ! for the patch with the most elements. - - if(allocated(site%taulambda_2str)) then - n_scr = ubound(site%taulambda_2str,dim=1) - allocate_scratch = .false. - if(2*max_elements > n_scr) then - allocate_scratch = .true. - deallocate(site%taulambda_2str,site%ipiv_2str,site%omega_2str) - elseif(2*max_elements < (n_scr-24)) then - allocate_scratch = .true. - deallocate(site%taulambda_2str,site%ipiv_2str,site%omega_2str) - end if - else - allocate_scratch = .true. - end if - - if(allocate_scratch)then - ! Twice as many spaces as there are elements, plus some - ! extra to prevent allocating/deallocating on the next step - n_scr = 2*max_elements+8 - allocate(site%taulambda_2str(n_scr)) - allocate(site%omega_2str(n_scr,n_scr)) - allocate(site%ipiv_2str(n_scr)) - end if - return end subroutine FatesConstructRadElements diff --git a/radiation/TwoStreamMLPEMod.F90 b/radiation/TwoStreamMLPEMod.F90 index 22bf303649..dfec3c1ebd 100644 --- a/radiation/TwoStreamMLPEMod.F90 +++ b/radiation/TwoStreamMLPEMod.F90 @@ -32,7 +32,7 @@ Module TwoStreamMLPEMod private real(r8),parameter :: nearzero = 1.e-20_r8 - logical, parameter :: debug=.true. + logical, parameter :: debug=.false. real(r8), parameter :: unset_r8 = 1.e-36_r8 real(r8), parameter :: unset_int = -999 integer, parameter :: twostr_vis = 1 ! Named index of visible shortwave radiation @@ -474,12 +474,20 @@ subroutine GetAbsRad(this,ican,icol,ib,vai_top,vai_bot, & real(r8) :: Rb_net ! Difference in beam radiation at upper and lower boundaries [W/m2] real(r8) :: vai_max ! total integrated (leaf+stem) area index of the current element real(r8) :: frac_abs_snow ! fraction of radiation absorbed by snow + real(r8) :: diff_wt_leaf ! diffuse absorption weighting for leaves real(r8) :: diff_wt_stem ! diffuse absorption weighting for stems + real(r8) :: diff_wt_snow ! diffuse absorption weighting for snow + real(r8) :: diff_wt_elem ! diffuse absorption normalization real(r8) :: beam_wt_leaf ! beam absorption weighting for leaves real(r8) :: beam_wt_stem ! beam absorption weighting for stems + real(r8) :: beam_wt_snow ! beam absorption weighting for snow + real(r8) :: beam_wt_elem ! beam absorption normalization + real(r8) :: lai_bot,lai_top real(r8) :: r_dn_top,r_dn_bot + + integer, parameter :: split_method = 2 associate(scelb => this%band(ib)%scelb(ican,icol), & scelg => this%scelg(ican,icol), & @@ -514,19 +522,6 @@ subroutine GetAbsRad(this,ican,icol,ib,vai_top,vai_bot, & leaf_sun_frac = 0001._r8 end if - ! We have to disentangle the absorption between leaves and stems, we give them both - ! a weighting fraction of total absorption of area*K*(1-om) - - frac_abs_snow = this%frac_snow*(1._r8-om_snow(ib)) / (1._r8-scelb%om) - - diff_wt_leaf = scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*rad_params%Kd_leaf(ft) - diff_wt_stem = scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*rad_params%Kd_stem(ft) - - beam_wt_leaf = scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*scelg%Kb_leaf - beam_wt_stem = scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*1._r8 - - ! Mean element transmission coefficients adding snow scattering - if(debug) then if( (vai_bot-vai_max)>rel_err_thresh)then write(log_unit,*)"During decomposition of the 2-stream radiation solution" @@ -546,7 +541,7 @@ subroutine GetAbsRad(this,ican,icol,ib,vai_top,vai_bot, & call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if - + ! Amount of absorbed radiation is retrieved by doing an energy ! balance on this boundaries over the depth of interest (ie net) ! Result is Watts / m2 of the element's area footprint NOT @@ -574,15 +569,80 @@ subroutine GetAbsRad(this,ican,icol,ib,vai_top,vai_bot, & Rb_abs = Rb_net * (1._r8-this%band(ib)%scelb(ican,icol)%om) Rd_abs = Rd_net + Rb_net * this%band(ib)%scelb(ican,icol)%om + + ! Disentangle the absorption between leaves and stems, we give them both + ! a weighting fraction of total absorption of area*K*(1-om) + + if(split_method == 1)then + + frac_abs_snow = this%frac_snow*(1._r8-om_snow(ib)) / (1._r8-scelb%om) + + diff_wt_leaf = scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*rad_params%Kd_leaf(ft) + diff_wt_stem = scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*rad_params%Kd_stem(ft) + + beam_wt_leaf = scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*scelg%Kb_leaf + beam_wt_stem = scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*1._r8 + Rb_abs_leaf = (1._r8-frac_abs_snow)*Rb_abs * beam_wt_leaf / (beam_wt_leaf+beam_wt_stem) + Rd_abs_leaf = (1._r8-frac_abs_snow)*Rd_abs * diff_wt_leaf / (diff_wt_leaf+diff_wt_stem) - Rb_abs_leaf = (1._r8-frac_abs_snow)*Rb_abs * beam_wt_leaf / (beam_wt_leaf+beam_wt_stem) - Rd_abs_leaf = (1._r8-frac_abs_snow)*Rd_abs * diff_wt_leaf / (diff_wt_leaf+diff_wt_stem) + R_abs_snow = (Rb_abs+Rd_abs)*frac_abs_snow + + R_abs_stem = (1._r8-frac_abs_snow)* & + (Rb_abs*beam_wt_stem / (beam_wt_leaf+beam_wt_stem) + & + Rd_abs*diff_wt_stem / (diff_wt_leaf+diff_wt_stem)) - R_abs_snow = (Rb_abs+Rd_abs)*frac_abs_snow + + elseif(split_method == 2)then - R_abs_stem = (1._r8-frac_abs_snow)* & - (Rb_abs*beam_wt_stem / (beam_wt_leaf+beam_wt_stem) + & - Rd_abs*diff_wt_stem / (diff_wt_leaf+diff_wt_stem)) + ! This is solely designed to have the same methods as how the cohort + ! caller disentangles absorptions and partitions into leaf,stem/snow + + diff_wt_leaf = (1._r8-this%frac_snow)*scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*rad_params%Kd_leaf(ft) + diff_wt_elem = (scelg%lai+scelg%sai)*(1._r8-scelb%om)*scelg%Kd + + beam_wt_leaf = (1._r8-this%frac_snow)*scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*scelg%Kb_leaf + beam_wt_elem = (scelg%lai+scelg%sai)*(1._r8-scelb%om)*scelg%Kb + + Rd_abs_leaf = Rd_abs * min(1.0_r8,diff_wt_leaf / diff_wt_elem) + Rb_abs_leaf = Rb_abs * min(1.0_r8,beam_wt_leaf / beam_wt_elem) + + R_abs_stem = (Rd_abs - Rd_abs_leaf)*(1._r8-this%frac_snow) + (Rb_abs-Rb_abs_leaf)*(1._r8-this%frac_snow) + + R_abs_snow = Rd_abs + Rb_abs - (Rd_abs_leaf+Rb_abs_leaf+R_abs_stem) + + else + + ! Improved Method + ! Calculate absorbance shares of Rb and Rd based off of their relative + ! optical depth (area index * K) times the material absorption (1-om) + ! Add up these shares to get the total for the element, use that as + ! a normalization factor + + diff_wt_snow = this%frac_snow*scelg%lai*(1._r8-om_snow(ib))*rad_params%Kd_leaf(ft) + diff_wt_snow = diff_wt_snow+this%frac_snow*scelg%sai*(1._r8-om_snow(ib))*rad_params%Kd_stem(ft) + diff_wt_leaf = (1._r8 - this%frac_snow)*scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*rad_params%Kd_leaf(ft) + diff_wt_stem = (1._r8 - this%frac_snow)*scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*rad_params%Kd_stem(ft) + diff_wt_elem = diff_wt_snow+diff_wt_leaf+diff_wt_stem + + + beam_wt_snow = this%frac_snow*scelg%lai*(1._r8-om_snow(ib))*scelg%Kb_leaf + beam_wt_snow = beam_wt_snow + this%frac_snow*scelg%sai*(1._r8-om_snow(ib))*1._r8 + beam_wt_leaf = (1._r8 - this%frac_snow)*scelg%lai*(1._r8-rad_params%om_leaf(ib,ft))*scelg%Kb_leaf + beam_wt_stem = (1._r8 - this%frac_snow)*scelg%sai*(1._r8-rad_params%om_stem(ib,ft))*1._r8 + beam_wt_elem = beam_wt_snow + beam_wt_leaf + beam_wt_stem + + + Rd_abs_leaf = Rd_abs * diff_wt_leaf/diff_wt_elem + Rb_abs_leaf = Rb_abs * beam_wt_leaf/beam_wt_elem + R_abs_stem = Rd_abs * diff_wt_stem/diff_wt_elem + Rb_abs * beam_wt_stem/beam_wt_elem + R_abs_snow = Rd_abs * diff_wt_snow/diff_wt_elem + Rb_abs * beam_wt_snow/beam_wt_elem + + + end if + + + + @@ -1103,9 +1163,6 @@ subroutine Solve(this, ib, & upper_boundary_type, & Rbeam_atm, & Rdiff_atm, & - taulamb, & - omega, & - ipiv, & albedo_beam, & albedo_diff, & consv_err, & @@ -1143,9 +1200,9 @@ subroutine Solve(this, ib, & real(r8) :: Rbeam_atm ! Intensity of beam radiation at top of canopy [W/m2 ground] real(r8) :: Rdiff_atm ! Intensity of diffuse radiation at top of canopy [W/m2 ground] ! - real(r8) :: taulamb(:) ! both the coefficient vector and constant side of the linear equation - real(r8) :: omega(:,:) ! the square matrix to be inverted - integer :: ipiv(:) ! pivot indices for LAPACK (not optional output, we don't use) + !real(r8) :: taulamb(:) ! both the coefficient vector and constant side of the linear equation + !real(r8) :: omega(:,:) ! the square matrix to be inverted + !integer :: ipiv(:) ! pivot indices for LAPACK (not optional output, we don't use) real(r8) :: albedo_beam ! Mean albedo at canopy top generated from beam radiation [-] real(r8) :: albedo_diff ! Mean albedo at canopy top generated from downwelling diffuse [-] @@ -1206,12 +1263,6 @@ subroutine Solve(this, ib, & ! Parameters for solving via LAPACK DGESV() and DGESVXX() integer :: info ! Procedure diagnostic ouput - ! Testing switch - ! If true, then allow elements - ! of different layers, but same row, to have priority - ! flux into the other element, instead of a mix - logical, parameter :: continuity_on = .true. - logical, parameter :: albedo_corr = .true. ! ------------------------------------------------------------------------------------ @@ -1341,11 +1392,13 @@ subroutine Solve(this, ib, & ! ===================================================================== n_eq = 2*this%n_scel - - ! TO-DO: MAKE THIS SCRATCH SPACE AT THE SITE LEVEL? - !!allocate(OMEGA(2*this%n_scel,2*this%n_scel),stat=alloc_err) - !!allocate(TAU(2*this%n_scel),stat=alloc_err) - !!allocate(LAMBDA(2*this%n_scel),stat=alloc_err) + + + block + + real(r8) :: omega(n_eq,n_eq) + real(r8) :: taulamb(n_eq) + integer :: ipiv(n_eq) ! We come up with two solutions: ! First: we run with now diffuse downwelling @@ -1573,22 +1626,23 @@ subroutine Solve(this, ib, & end if ! Perform a forward check on the solution error - do ilem = 1,n_eq - temp_err = tau_temp(ilem) - sum(taulamb(1:n_eq)*omega_temp(ilem,1:n_eq)) - if(abs(temp_err)>rel_err_thresh)then - write(log_unit,*) 'Poor forward solution on two-stream solver' - write(log_unit,*) 'isol (1=beam or 2=diff): ',isol - write(log_unit,*) 'i (equation): ',ilem - write(log_unit,*) 'band index (1=vis,2=nir): ',ib - write(log_unit,*) 'error (tau(i) - omega(i,:)*lambda(:)) ',temp_err - this%band(ib)%Rbeam_atm = 1._r8 - this%band(ib)%Rdiff_atm = 1._r8 - call this%Dump(ib) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - end do - deallocate(tau_temp,omega_temp) - + if(debug)then + do ilem = 1,n_eq + temp_err = tau_temp(ilem) - sum(taulamb(1:n_eq)*omega_temp(ilem,1:n_eq)) + if(abs(temp_err)>rel_err_thresh)then + write(log_unit,*) 'Poor forward solution on two-stream solver' + write(log_unit,*) 'isol (1=beam or 2=diff): ',isol + write(log_unit,*) 'i (equation): ',ilem + write(log_unit,*) 'band index (1=vis,2=nir): ',ib + write(log_unit,*) 'error (tau(i) - omega(i,:)*lambda(:)) ',temp_err + this%band(ib)%Rbeam_atm = 1._r8 + this%band(ib)%Rdiff_atm = 1._r8 + call this%Dump(ib) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + deallocate(tau_temp,omega_temp) + end if ! Save the solution terms @@ -1702,6 +1756,8 @@ subroutine Solve(this, ib, & end do do_isol + end block + ! Check the error balance ! ---------------------------------------------------------------------------------------------