diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 90b7136444..f5ec118c24 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1229,7 +1229,7 @@ subroutine canopy_spread( currentSite ) do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) - if( ( int(prt_params%woody(currentCohort%pft)) .eq. itrue ) .and. & + if( (prt_params%woody(currentCohort%pft) .eq. itrue ) .and. & (currentCohort%canopy_layer .eq. 1 ) ) then sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area endif @@ -1340,7 +1340,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) if(currentCohort%canopy_layer==1)then currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area - if( int(prt_params%woody(ft))==itrue)then + if( prt_params%woody(ft)==itrue)then currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area endif endif diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 4318ee469f..7e10496679 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -70,6 +70,7 @@ module EDCohortDynamicsMod use FatesAllometryMod , only : tree_lai, tree_sai use FatesAllometryMod , only : set_root_fraction use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes use PRTGenericMod, only : all_carbon_elements @@ -85,6 +86,7 @@ module EDCohortDynamicsMod use PRTGenericMod, only : SetState use PRTAllometricCarbonMod, only : callom_prt_vartypes + use PRTAllometricCarbonMod, only : csimpler_allom_prt_vartypes use PRTAllometricCarbonMod, only : ac_bc_inout_id_netdc use PRTAllometricCarbonMod, only : ac_bc_in_id_pft use PRTAllometricCarbonMod, only : ac_bc_in_id_ctrim @@ -147,7 +149,7 @@ module EDCohortDynamicsMod subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & - prt, laimemory, sapwmemory, structmemory, & + prt, leafmemory, sapwmemory, structmemory, & status, recruitstatus,ctrim, carea, clayer, spread, bc_in) ! ! !DESCRIPTION: @@ -181,7 +183,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & real(r8), intent(in) :: dbh ! dbh: cm class(prt_vartypes),target :: prt ! The allocated PARTEH ! object - real(r8), intent(in) :: laimemory ! target leaf biomass- set from + real(r8), intent(in) :: leafmemory ! target leaf biomass- set from ! previous year: kGC per indiv real(r8), intent(in) :: sapwmemory ! target sapwood biomass- set from ! previous year: kGC per indiv @@ -237,7 +239,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & new_cohort%canopy_trim = ctrim new_cohort%canopy_layer = clayer new_cohort%canopy_layer_yesterday = real(clayer, r8) - new_cohort%laimemory = laimemory + new_cohort%leafmemory = leafmemory new_cohort%sapwmemory = sapwmemory new_cohort%structmemory = structmemory @@ -400,7 +402,7 @@ subroutine InitPRTBoundaryConditions(new_cohort) select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) ! Register boundary conditions for the Carbon Only Allometric Hypothesis @@ -463,10 +465,16 @@ subroutine InitPRTObject(prt) ! Potential Extended types class(callom_prt_vartypes), pointer :: c_allom_prt + class(csimpler_allom_prt_vartypes), pointer :: csimpler_allom_prt class(cnp_allom_prt_vartypes), pointer :: cnp_allom_prt select case(hlm_parteh_mode) + case (prt_csimpler_allom_hyp) + + allocate(csimpler_allom_prt) + prt => csimpler_allom_prt + case (prt_carbon_allom_hyp) allocate(c_allom_prt) @@ -540,8 +548,8 @@ subroutine nan_cohort(cc_p) currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm currentCohort%coage = nan ! age of the cohort in years - currentCohort%hite = nan ! height: meters - currentCohort%laimemory = nan ! target leaf biomass- set from previous year: kGC per indiv + currentCohort%hite = nan ! height: meters + currentCohort%leafmemory = nan ! target leaf biomass- set from previous year: kGC per indiv currentCohort%sapwmemory = nan ! target sapwood biomass- set from previous year: kGC per indiv currentCohort%structmemory = nan ! target structural biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 @@ -765,7 +773,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ if (currentcohort%n < min_n_safemath .and. level == 1) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,call_index + write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index endif endif @@ -778,7 +786,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ (currentCohort%dbh < 0.00001_r8 .and. store_c < 0._r8) ) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,call_index + write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index endif endif @@ -786,7 +794,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ if (currentCohort%canopy_layer > nclmax ) then terminate = itrue if ( debug ) then - write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,call_index + write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index endif endif @@ -796,7 +804,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ terminate = itrue if ( debug ) then write(fates_log(),*) 'terminating cohorts 3', & - sapw_c,leaf_c,fnrt_c,store_c,call_index + sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index endif endif @@ -805,7 +813,7 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_ terminate = itrue if ( debug ) then write(fates_log(),*) 'terminating cohorts 4', & - struct_c,sapw_c,leaf_c,fnrt_c,store_c,call_index + struct_c,sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index endif endif @@ -1196,7 +1204,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'Cohort I, Cohort II' write(fates_log(),*) 'n:',currentCohort%n,nextc%n write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew - write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory + write(fates_log(),*) 'leafmemory:',currentCohort%leafmemory,nextc%leafmemory write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh @@ -1234,8 +1242,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ----------------------------------------------------------------- call UpdateCohortBioPhysRates(currentCohort) - currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & - + nextc%n*nextc%laimemory)/newn + currentCohort%leafmemory = (currentCohort%n*currentCohort%leafmemory & + + nextc%n*nextc%leafmemory)/newn currentCohort%sapwmemory = (currentCohort%n*currentCohort%sapwmemory & + nextc%n*nextc%sapwmemory)/newn @@ -1828,7 +1836,7 @@ subroutine copy_cohort( currentCohort,copyc ) n%dbh = o%dbh n%coage = o%coage n%hite = o%hite - n%laimemory = o%laimemory + n%leafmemory = o%leafmemory n%sapwmemory = o%sapwmemory n%structmemory = o%structmemory n%lai = o%lai @@ -2107,7 +2115,7 @@ subroutine EvaluateAndCorrectDBH(currentCohort,delta_dbh,delta_hite) delta_dbh = 0._r8 delta_hite = 0._r8 - if( int(prt_params%woody(currentCohort%pft)) == itrue) then + if( prt_params%woody(currentCohort%pft) == itrue) then struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index f1f23d9f33..1da4ae4d6b 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -258,7 +258,7 @@ subroutine LoggingMortality_frac( pft_i, dbh, canopy_layer, lmort_direct, & ! transfer of area to secondary land is based on overall area affected, not just logged crown area ! l_degrad accounts for the affected area between logged crowns - if(int(prt_params%woody(pft_i)) == 1)then ! only set logging rates for trees + if(prt_params%woody(pft_i) == 1)then ! only set logging rates for trees ! direct logging rates, based on dbh min and max criteria if (dbh >= logging_dbhmin .and. .not. & @@ -542,7 +542,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site ! plants that were impacted. Thus, no direct dead can occur ! here, and indirect are impacts. - if(int(prt_params%woody(pft)) == itrue) then + if(prt_params%woody(pft) == itrue) then direct_dead = 0.0_r8 indirect_dead = logging_coll_under_frac * & (1._r8-currentPatch%fract_ldist_not_harvested) * currentCohort%n * & diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 537e74824d..a6c4ef7d02 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -84,7 +84,8 @@ module EDPatchDynamicsMod use PRTGenericMod, only : struct_organ use PRTLossFluxesMod, only : PRTBurnLosses use FatesInterfaceTypesMod, only : hlm_parteh_mode - use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use SFParamsMod, only : SF_VAL_CWD_FRAC use EDParamsMod, only : logging_event_code @@ -777,7 +778,7 @@ subroutine spawn_patches( currentSite, bc_in) else ! small trees - if( int(prt_params%woody(currentCohort%pft)) == itrue)then + if( prt_params%woody(currentCohort%pft) == itrue)then ! Survivorship of undestory woody plants. Two step process. @@ -917,7 +918,7 @@ subroutine spawn_patches( currentSite, bc_in) ! burned off. Here, we remove that mass, and ! tally it in the flux we sent to the atmosphere - if(int(prt_params%woody(currentCohort%pft)) == itrue)then + if( prt_params%woody(currentCohort%pft) == itrue)then leaf_burn_frac = currentCohort%fraction_crown_burned else @@ -995,7 +996,7 @@ subroutine spawn_patches( currentSite, bc_in) ! WHat to do with cohorts in the understory of a logging generated ! disturbance patch? - if(int(prt_params%woody(currentCohort%pft)) == itrue)then + if(prt_params%woody(currentCohort%pft) == itrue)then ! Survivorship of undestory woody plants. Two step process. @@ -1905,7 +1906,7 @@ subroutine mortality_litter_fluxes(currentSite, currentPatch, & num_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * & hlm_freq_day * fates_mortality_disturbance_fraction) - elseif(int(prt_params%woody(pft)) == itrue) then + elseif(prt_params%woody(pft) == itrue) then ! Understorey trees. The total dead is based on their survivorship ! function, and the total area of disturbance. diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 5d66f56d39..c4e8085fba 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -87,6 +87,7 @@ module EDPhysiologyMod use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : CheckIntegratedAllometries use FatesAllometryMod, only : set_root_fraction + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -424,9 +425,9 @@ subroutine trim_canopy( currentSite ) real(r8) :: work(workmax) ! work array real(r8) :: initial_trim ! Initial trim - real(r8) :: optimum_trim ! Optimum trim value - real(r8) :: initial_laimem ! Initial laimemory - real(r8) :: optimum_laimem ! Optimum laimemory + real(r8) :: optimum_trim ! Optimum trim value + real(r8) :: initial_leafmem ! Initial leafmemory + real(r8) :: optimum_leafmem ! Optimum leafmemory !---------------------------------------------------------------------- @@ -446,15 +447,15 @@ subroutine trim_canopy( currentSite ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) - ! Save off the incoming trim and laimemory + ! Save off the incoming trim and leafmemory initial_trim = currentCohort%canopy_trim - initial_laimem = currentCohort%laimemory + initial_leafmem = currentCohort%leafmemory ! Add debug diagnstic output to determine which cohort if (debug) then write(fates_log(),*) 'Current cohort:', icohort write(fates_log(),*) 'Starting canopy trim:', initial_trim - write(fates_log(),*) 'Starting laimemory:', currentCohort%laimemory + write(fates_log(),*) 'Starting leafmemory:', currentCohort%leafmemory endif trimmed = .false. @@ -601,7 +602,7 @@ subroutine trim_canopy( currentSite ) currentCohort%canopy_trim = currentCohort%canopy_trim - & EDPftvarcon_inst%trim_inc(ipft) if (prt_params%evergreen(ipft) /= 1)then - currentCohort%laimemory = currentCohort%laimemory * & + currentCohort%leafmemory = currentCohort%leafmemory * & (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) endif @@ -649,15 +650,15 @@ subroutine trim_canopy( currentSite ) ! optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim - optimum_laimem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_laimem + optimum_leafmem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_leafmem ! Determine if the optimum trim value makes sense. The smallest cohorts tend to have unrealistic fits. if (optimum_trim > 0. .and. optimum_trim < 1.) then currentCohort%canopy_trim = optimum_trim - ! If the cohort pft is not evergreen we reduce the laimemory as well + ! If the cohort pft is not evergreen we reduce the leafmemory as well if (prt_params%evergreen(ipft) /= 1) then - currentCohort%laimemory = optimum_laimem + currentCohort%leafmemory = optimum_leafmem endif trimmed = .true. @@ -1123,14 +1124,14 @@ subroutine phenology_leafonoff(currentSite) ! stop flow of carbon out of bstore. if(store_c>nearzero) then - ! flush either the amount required from the laimemory, or -most- of the storage pool + ! flush either the amount required from the leafmemory, or -most- of the storage pool ! RF: added a criterion to stop the entire store pool emptying and triggering termination mortality ! n.b. this might not be necessary if we adopted a more gradual approach to leaf flushing... store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & - currentCohort%laimemory)/store_c,(1.0_r8-carbon_store_buffer)) + currentCohort%leafmemory)/store_c,(1.0_r8-carbon_store_buffer)) if(prt_params%woody(ipft).ne.itrue)then - totalmemory=currentCohort%laimemory+currentCohort%sapwmemory+currentCohort%structmemory + totalmemory=currentCohort%leafmemory+currentCohort%sapwmemory+currentCohort%structmemory store_c_transfer_frac = min((EDPftvarcon_inst%phenflush_fraction(ipft)* & totalmemory)/store_c, (1.0_r8-carbon_store_buffer)) endif @@ -1144,7 +1145,7 @@ subroutine phenology_leafonoff(currentSite) if(prt_params%woody(ipft) == itrue) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, store_c_transfer_frac) - currentCohort%laimemory = 0.0_r8 + currentCohort%leafmemory = 0.0_r8 else @@ -1152,7 +1153,7 @@ subroutine phenology_leafonoff(currentSite) if (stem_drop_fraction .gt. 0.0_r8) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%laimemory/totalmemory) + store_c_transfer_frac*currentCohort%leafmemory/totalmemory) call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) @@ -1167,7 +1168,7 @@ subroutine phenology_leafonoff(currentSite) end if - currentCohort%laimemory = 0.0_r8 + currentCohort%leafmemory = 0.0_r8 currentCohort%structmemory = 0.0_r8 currentCohort%sapwmemory = 0.0_r8 @@ -1188,10 +1189,9 @@ subroutine phenology_leafonoff(currentSite) ! This sets the cohort to the "leaves off" flag currentCohort%status_coh = leaves_off - ! Remember what the lai was (leaf mass actually) was for next year + ! Remember what the leaf mass was for next year ! the same amount back on in the spring... - - currentCohort%laimemory = leaf_c + currentCohort%leafmemory = leaf_c ! Drop Leaves (this routine will update the leaf state variables, ! for carbon and any other element that are prognostic. It will @@ -1238,12 +1238,12 @@ subroutine phenology_leafonoff(currentSite) if(store_c>nearzero) then store_c_transfer_frac = & - min((EDPftvarcon_inst%phenflush_fraction(ipft)*currentCohort%laimemory)/store_c, & + min((EDPftvarcon_inst%phenflush_fraction(ipft)*currentCohort%leafmemory)/store_c, & (1.0_r8-carbon_store_buffer)) if(prt_params%woody(ipft).ne.itrue)then - totalmemory=currentCohort%laimemory+currentCohort%sapwmemory+currentCohort%structmemory + totalmemory=currentCohort%leafmemory+currentCohort%sapwmemory+currentCohort%structmemory store_c_transfer_frac = min(EDPftvarcon_inst%phenflush_fraction(ipft)*totalmemory/store_c, & (1.0_r8-carbon_store_buffer)) @@ -1260,7 +1260,7 @@ subroutine phenology_leafonoff(currentSite) call PRTPhenologyFlush(currentCohort%prt, ipft, & leaf_organ, store_c_transfer_frac) - currentCohort%laimemory = 0.0_r8 + currentCohort%leafmemory = 0.0_r8 else @@ -1268,7 +1268,7 @@ subroutine phenology_leafonoff(currentSite) if (stem_drop_fraction .gt. 0.0_r8) then call PRTPhenologyFlush(currentCohort%prt, ipft, leaf_organ, & - store_c_transfer_frac*currentCohort%laimemory/totalmemory) + store_c_transfer_frac*currentCohort%leafmemory/totalmemory) call PRTPhenologyFlush(currentCohort%prt, ipft, sapw_organ, & store_c_transfer_frac*currentCohort%sapwmemory/totalmemory) @@ -1283,7 +1283,7 @@ subroutine phenology_leafonoff(currentSite) end if - currentCohort%laimemory = 0.0_r8 + currentCohort%leafmemory = 0.0_r8 currentCohort%structmemory = 0.0_r8 currentCohort%sapwmemory = 0.0_r8 @@ -1300,8 +1300,8 @@ subroutine phenology_leafonoff(currentSite) ! This sets the cohort to the "leaves off" flag currentCohort%status_coh = leaves_off - ! Remember what the lai (leaf mass actually) was for next year - currentCohort%laimemory = leaf_c + ! Remember what the leaf mass was for next year + currentCohort%leafmemory = leaf_c call PRTDeciduousTurnover(currentCohort%prt,ipft, & leaf_organ, leaf_drop_fraction) @@ -1890,7 +1890,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! Default assumption is that leaves are on cohortstatus = leaves_on - temp_cohort%laimemory = 0.0_r8 + temp_cohort%leafmemory = 0.0_r8 temp_cohort%sapwmemory = 0.0_r8 temp_cohort%structmemory = 0.0_r8 @@ -1899,7 +1899,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! as "cold", then set the cohort's status to leaves_off, and remember the leaf biomass if ((prt_params%season_decid(ft) == itrue) .and. & (any(currentSite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold]))) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf c_leaf = 0.0_r8 ! If plant is not woody then set sapwood and structural biomass as well @@ -1917,7 +1917,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! biomass if ((prt_params%stress_decid(ft) == itrue) .and. & (any(currentSite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff]))) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf c_leaf = 0.0_r8 ! If plant is not woody then set sapwood and structural biomass as well @@ -2044,7 +2044,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Put all of the leaf mass into the first bin call SetState(prt,leaf_organ, element_id,m_leaf,1) @@ -2103,7 +2103,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! This initializes the cohort call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, & temp_cohort%hite, temp_cohort%coage, temp_cohort%dbh, prt, & - temp_cohort%laimemory, temp_cohort%sapwmemory, temp_cohort%structmemory, & + temp_cohort%leafmemory, temp_cohort%sapwmemory, temp_cohort%structmemory, & cohortstatus, recruitstatus, & temp_cohort%canopy_trim,temp_cohort%c_area, & currentPatch%NCL_p, currentSite%spread, bc_in) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 42264ca776..dc767d08f4 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2370,7 +2370,7 @@ subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bl ) integer, parameter :: max_counter = 200 ! Do reduce "if" calls, we break this call into two parts - if ( int(prt_params%woody(ipft)) == itrue ) then + if ( prt_params%woody(ipft) == itrue ) then if(.not.present(bdead)) then write(fates_log(),*) 'woody plants must use structure for dbh reset' @@ -2456,7 +2456,7 @@ subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bl ) call h_allom(d,ipft,h) if(counter>10)then write(fates_log(),*) 'dbh counter: ',counter,' is woody: ',& - int(prt_params%woody(ipft))==itrue + prt_params%woody(ipft)==itrue end if diff --git a/biogeochem/FatesSoilBGCFluxMod.F90 b/biogeochem/FatesSoilBGCFluxMod.F90 index d14ce7b005..c6d7089802 100644 --- a/biogeochem/FatesSoilBGCFluxMod.F90 +++ b/biogeochem/FatesSoilBGCFluxMod.F90 @@ -15,6 +15,7 @@ module FatesSoilBGCFluxMod use PRTGenericMod , only : num_elements use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : prt_vartypes @@ -219,7 +220,8 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) end do ! We can exit if this is a c-only simulation - if(hlm_parteh_mode.eq.prt_carbon_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) ! These can now be zero'd do s = 1, nsites bc_in(s)%plant_nh4_uptake_flux(:,:) = 0._r8 @@ -227,7 +229,7 @@ subroutine UnPackNutrientAquisitionBCs(sites, bc_in) bc_in(s)%plant_p_uptake_flux(:,:) = 0._r8 end do return - end if + end select do s = 1, nsites diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90 index adffe67d85..54138ffaaf 100644 --- a/biogeophys/FatesPlantRespPhotosynthMod.F90 +++ b/biogeophys/FatesPlantRespPhotosynthMod.F90 @@ -39,6 +39,7 @@ module FATESPlantRespPhotosynthMod use PRTGenericMod, only : max_nleafage use EDTypesMod, only : do_fates_salinity use EDParamsMod, only : q10_mr + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : all_carbon_elements @@ -417,7 +418,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) if ( .not.rate_mask_z(iv,ft,cl) .or. & (hlm_use_planthydro.eq.itrue) .or. & (nleafage > 1) .or. & - (hlm_parteh_mode .ne. prt_carbon_allom_hyp ) ) then + (hlm_parteh_mode .ne. prt_csimpler_allom_hyp .and. & + hlm_parteh_mode .ne. prt_carbon_allom_hyp) ) then if (hlm_use_planthydro.eq.itrue ) then @@ -480,7 +482,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) ! 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) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) lnc_top = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft) @@ -500,6 +502,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) end select + ! MLO - Shouldn't these numbers be parameters too? lmr25top = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8)) lmr25top = lmr25top * lnc_top / (umolC_to_kgC * g_per_kg) @@ -639,7 +642,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime) fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) 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)) diff --git a/fire/SFMainMod.F90 b/fire/SFMainMod.F90 index aedcb4aa7c..ee8e4c2400 100644 --- a/fire/SFMainMod.F90 +++ b/fire/SFMainMod.F90 @@ -192,7 +192,7 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%livegrass = 0.0_r8 currentCohort => currentPatch%tallest do while(associated(currentCohort)) - if( int(prt_params%woody(currentCohort%pft)) == ifalse)then + if( prt_params%woody(currentCohort%pft) == ifalse)then currentPatch%livegrass = currentPatch%livegrass + & currentCohort%prt%GetState(leaf_organ, all_carbon_elements) * & @@ -374,7 +374,7 @@ subroutine wind_effect ( currentSite, bc_in) do while(associated(currentCohort)) if (debug) write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area - if( int(prt_params%woody(currentCohort%pft)) == itrue)then + if( prt_params%woody(currentCohort%pft) == itrue)then currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area else total_grass_area = total_grass_area + currentCohort%c_area @@ -864,7 +864,7 @@ subroutine crown_scorching ( currentSite ) if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; do while(associated(currentCohort)) - if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only + if ( prt_params%woody(currentCohort%pft) == itrue) then !trees only leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) @@ -878,7 +878,7 @@ subroutine crown_scorching ( currentSite ) enddo !end cohort loop do i_pft=1,numpft - if (tree_ag_biomass > 0.0_r8 .and. int(prt_params%woody(i_pft)) == itrue) then + if (tree_ag_biomass > 0.0_r8 .and. prt_params%woody(i_pft) == itrue) then !Equation 16 in Thonicke et al. 2010 !Van Wagner 1973 EQ8 !2/3 Byram (1959) currentPatch%Scorch_ht(i_pft) = EDPftvarcon_inst%fire_alpha_SH(i_pft) * (currentPatch%FI**0.667_r8) @@ -920,7 +920,7 @@ subroutine crown_damage ( currentSite ) do while(associated(currentCohort)) currentCohort%fraction_crown_burned = 0.0_r8 - if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only + if ( prt_params%woody(currentCohort%pft) == itrue) then !trees only ! Flames lower than bottom of canopy. ! c%hite is height of cohort @@ -984,7 +984,7 @@ subroutine cambial_damage_kill ( currentSite ) if (currentPatch%fire == 1) then currentCohort => currentPatch%tallest; do while(associated(currentCohort)) - if ( int(prt_params%woody(currentCohort%pft)) == itrue) then !trees only + if ( prt_params%woody(currentCohort%pft) == itrue) then !trees only ! Equation 21 in Thonicke et al 2010 bt = EDPftvarcon_inst%bark_scaler(currentCohort%pft)*currentCohort%dbh ! bark thickness. ! Equation 20 in Thonicke et al. 2010. @@ -1036,7 +1036,7 @@ subroutine post_fire_mortality ( currentSite ) do while(associated(currentCohort)) currentCohort%fire_mort = 0.0_r8 currentCohort%crownfire_mort = 0.0_r8 - if ( int(prt_params%woody(currentCohort%pft)) == itrue) then + if ( prt_params%woody(currentCohort%pft) == itrue) then ! Equation 22 in Thonicke et al. 2010. currentCohort%crownfire_mort = EDPftvarcon_inst%crown_kill(currentCohort%pft)*currentCohort%fraction_crown_burned**3.0_r8 ! Equation 18 in Thonicke et al. 2010. diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 06bcd1858a..d894978f3a 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -61,6 +61,7 @@ module EDInitMod use FatesAllometryMod , only : bstore_allom use PRTGenericMod , only : StorageNutrientTarget use FatesInterfaceTypesMod, only : hlm_parteh_mode + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -787,7 +788,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call bstore_allom(temp_cohort%dbh, pft, temp_cohort%canopy_trim, c_store) - temp_cohort%laimemory = 0._r8 + temp_cohort%leafmemory = 0._r8 temp_cohort%sapwmemory = 0._r8 temp_cohort%structmemory = 0._r8 cstatus = leaves_on @@ -798,7 +799,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) if( prt_params%season_decid(pft) == itrue .and. & any(site_in%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf temp_cohort%sapwmemory = c_sapw * stem_drop_fraction temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 @@ -809,7 +810,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) if ( prt_params%stress_decid(pft) == itrue .and. & any(site_in%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf temp_cohort%sapwmemory = c_sapw * stem_drop_fraction temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 @@ -850,26 +851,26 @@ subroutine init_cohorts( site_in, patch_in, bc_in) case(nitrogen_element) - m_struct = c_struct*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) - m_leaf = c_leaf*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) - m_fnrt = c_fnrt*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) - m_sapw = c_sapw*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) + m_struct = c_struct*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) + m_leaf = c_leaf*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) + m_fnrt = c_fnrt*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) + m_sapw = c_sapw*prt_params%nitr_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) m_repro = 0._r8 - m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) + m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) case(phosphorus_element) - m_struct = c_struct*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) - m_leaf = c_leaf*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) - m_fnrt = c_fnrt*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) - m_sapw = c_sapw*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) + m_struct = c_struct*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(struct_organ)) + m_leaf = c_leaf*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(leaf_organ)) + m_fnrt = c_fnrt*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(fnrt_organ)) + m_sapw = c_sapw*prt_params%phos_stoich_p2(pft,prt_params%organ_param_id(sapw_organ)) m_repro = 0._r8 - m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) + m_store = StorageNutrientTarget(pft,element_id,m_leaf,m_fnrt,m_sapw,m_struct) end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Put all of the leaf mass into the first bin call SetState(prt_obj,leaf_organ, element_id,m_leaf,1) @@ -893,7 +894,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call prt_obj%CheckInitialConditions() call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, & - temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%laimemory, & + temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%leafmemory, & temp_cohort%sapwmemory, temp_cohort%structmemory, cstatus, rstatus, & temp_cohort%canopy_trim, temp_cohort%c_area,1, site_in%spread, bc_in) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index a149132a8d..a6c70ee18f 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -17,7 +17,9 @@ module EDPftvarcon use FatesLitterMod, only : ilabile,icellulose,ilignin use PRTGenericMod, only : leaf_organ, fnrt_organ, store_organ use PRTGenericMod, only : sapw_organ, struct_organ, repro_organ - use PRTGenericMod, only : prt_cnp_flex_allom_hyp,prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp use FatesInterfaceTypesMod, only : hlm_nitrogen_spec, hlm_phosphorus_spec use FatesInterfaceTypesMod, only : hlm_parteh_mode use FatesInterfaceTypesMod, only : hlm_nu_com @@ -1502,7 +1504,8 @@ subroutine FatesCheckParams(is_master) if(.not.is_master) return - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) ! Check to see if either RD/ECA/MIC is turned on @@ -1537,15 +1540,19 @@ subroutine FatesCheckParams(is_master) end if end if + + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp) + ! No additional checks needed for now. + continue - elseif (hlm_parteh_mode .ne. prt_carbon_allom_hyp) then + case default write(fates_log(),*) 'FATES Plant Allocation and Reactive Transport has' write(fates_log(),*) 'only 2 modules supported, allometric carbon and CNP.' write(fates_log(),*) 'fates_parteh_mode must be set to 1 or 2 in the namelist' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + end select ! If any PFTs are specified as either prescribed N or P uptake ! then they all must be ! diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9ce4d5fe44..d377002aa3 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -219,7 +219,7 @@ module EDTypesMod real(r8) :: coage ! cohort age in years real(r8) :: hite ! height: meters integer :: indexnumber ! unique number for each cohort. (within clump?) - real(r8) :: laimemory ! target leaf biomass- set from previous year: kGC per indiv + real(r8) :: leafmemory ! target leaf biomass- set from previous year: kGC per indiv real(r8) :: sapwmemory ! target sapwood biomass- set from previous year: kGC per indiv real(r8) :: structmemory ! target structural biomass- set from previous year: kGC per indiv integer :: canopy_layer ! canopy status of cohort (1 = canopy, 2 = understorey, etc.) @@ -1049,7 +1049,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%dbh = ', ccohort%dbh write(fates_log(),*) 'co%hite = ', ccohort%hite write(fates_log(),*) 'co%coage = ', ccohort%coage - write(fates_log(),*) 'co%laimemory = ', ccohort%laimemory + write(fates_log(),*) 'co%leafmemory = ', ccohort%leafmemory write(fates_log(),*) 'co%sapwmemory = ', ccohort%sapwmemory write(fates_log(),*) 'co%structmemory = ', ccohort%structmemory diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index f43289da15..784892fea4 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -78,6 +78,7 @@ module FatesHistoryInterfaceMod use PRTGenericMod , only : all_carbon_elements use PRTGenericMod , only : carbon12_element use PRTGenericMod , only : nitrogen_element, phosphorus_element + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp implicit none @@ -221,6 +222,9 @@ module FatesHistoryInterfaceMod integer :: ih_bleaf_canopy_si_scpf integer :: ih_bleaf_understory_si_scpf + ! Size-class x PFT LAI states + integer :: ih_lai_canopy_si_scpf + integer :: ih_lai_understory_si_scpf integer :: ih_totvegn_scpf @@ -1861,6 +1865,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_bstor_understory_si_scpf => this%hvars(ih_bstor_understory_si_scpf)%r82d, & hio_bleaf_canopy_si_scpf => this%hvars(ih_bleaf_canopy_si_scpf)%r82d, & hio_bleaf_understory_si_scpf => this%hvars(ih_bleaf_understory_si_scpf)%r82d, & + hio_lai_canopy_si_scpf => this%hvars(ih_lai_canopy_si_scpf)%r82d, & + hio_lai_understory_si_scpf => this%hvars(ih_lai_understory_si_scpf)%r82d, & hio_mortality_canopy_si_scpf => this%hvars(ih_mortality_canopy_si_scpf)%r82d, & hio_mortality_understory_si_scpf => this%hvars(ih_mortality_understory_si_scpf)%r82d, & hio_nplant_canopy_si_scpf => this%hvars(ih_nplant_canopy_si_scpf)%r82d, & @@ -2506,7 +2512,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) store_m_net_alloc*n_perm2 / days_per_year / sec_per_day ! Woody State Variables (basal area growth increment) - if ( int(prt_params%woody(ft)) == itrue) then + if ( prt_params%woody(ft) == itrue) then ! basal area [m2/m2] hio_ba_si_scpf(io_si,scpf) = hio_ba_si_scpf(io_si,scpf) + & @@ -2633,6 +2639,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) store_m * ccohort%n / m2_per_ha hio_bleaf_canopy_si_scpf(io_si,scpf) = hio_bleaf_canopy_si_scpf(io_si,scpf) + & leaf_m * ccohort%n / m2_per_ha + hio_lai_canopy_si_scpf(io_si,scpf) = hio_lai_canopy_si_scpf(io_si,scpf) + & + ccohort%treelai*ccohort%c_area * AREA_INV hio_canopy_biomass_si(io_si) = hio_canopy_biomass_si(io_si) + n_perm2 * total_m @@ -2726,6 +2734,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) leaf_m * ccohort%n / m2_per_ha hio_understory_biomass_si(io_si) = hio_understory_biomass_si(io_si) + & n_perm2 * total_m + hio_lai_understory_si_scpf(io_si,scpf) = hio_lai_understory_si_scpf(io_si,scpf) + & + ccohort%treelai*ccohort%c_area * AREA_INV !hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + @@ -4044,7 +4054,8 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,bc_in,dt_tstep) do j_bc = j_t,j_b vwc = bc_in(s)%h2o_liqvol_sl(j_bc) - psi = site_hydr%wrf_soil(j)%p%psi_from_th(vwc) + psi = site_hydr%wrf_soil(j)%p%psi_from_th(vwc) ! MLO: Any reason for not using smp_sl? + ! cap capillary pressure ! psi = max(-1e5_r8,psi) Removing cap as that is inconstistent ! with model internals and physics. Should @@ -5792,6 +5803,14 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', upfreq=1, & ivar=ivar, initialize=initialize_variables, & index = ih_bleaf_canopy_si_scpf) + + call this%set_history_var(vname='FATES_LAI_CANOPY_SZPF', & + units = 'm2 m-2', & + long='Leaf area index (LAI) of canopy plants by pft/size', & + use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', upfreq=1, & + ivar=ivar, initialize=initialize_variables, & + index = ih_lai_canopy_si_scpf ) call this%set_history_var(vname='FATES_NPLANT_CANOPY_SZPF', units = 'm-2', & long='number of canopy plants by size/pft per m2', & @@ -5820,6 +5839,13 @@ subroutine define_history_vars(this, initialize_variables) use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & hlms='CLM:ALM', upfreq=1, ivar=ivar, & initialize=initialize_variables, index = ih_bleaf_understory_si_scpf) + + call this%set_history_var(vname='FATES_LAI_UNDERSTORY_SZPF', & + units = 'm2 m-2', & + long='Leaf area index (LAI) of understory plants by pft/size', & + use_default='inactive', avgflag='A', vtype=site_size_pft_r8, & + hlms='CLM:ALM', upfreq=1, ivar=ivar, & + initialize=initialize_variables, index = ih_lai_understory_si_scpf ) call this%set_history_var(vname='FATES_NPLANT_USTORY_SZPF', & units = 'm-2', & diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 3ce7589f48..d9be9c43f6 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -65,6 +65,7 @@ module FatesInterfaceMod use PRTGenericMod , only : element_list use PRTGenericMod , only : element_pos use EDParamsMod , only : eca_plant_escalar + use PRTGenericMod , only : prt_csimpler_allom_hyp use PRTGenericMod , only : prt_carbon_allom_hyp use PRTGenericMod , only : prt_cnp_flex_allom_hyp use PRTGenericMod , only : carbon12_element @@ -306,7 +307,7 @@ subroutine zero_bcs(fates,s) ! Fates -> BGC fragmentation mass fluxes select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) fates%bc_out(s)%litt_flux_cel_c_si(:) = 0._r8 fates%bc_out(s)%litt_flux_lig_c_si(:) = 0._r8 fates%bc_out(s)%litt_flux_lab_c_si(:) = 0._r8 @@ -627,7 +628,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in) ! Fates -> BGC fragmentation mass fluxes select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) allocate(bc_out%litt_flux_cel_c_si(nlevdecomp_in)) allocate(bc_out%litt_flux_lig_c_si(nlevdecomp_in)) allocate(bc_out%litt_flux_lab_c_si(nlevdecomp_in)) @@ -939,7 +940,7 @@ subroutine InitPARTEHGlobals() ! automatically. select case(hlm_parteh_mode) - case(prt_carbon_allom_hyp) + case(prt_csimpler_allom_hyp,prt_carbon_allom_hyp) num_elements = 1 allocate(element_list(num_elements)) diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 507f01dbee..a8e744446e 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -48,6 +48,7 @@ module FatesInventoryInitMod use EDPftvarcon , only : EDPftvarcon_inst use FatesInterfaceTypesMod, only : hlm_parteh_mode use EDCohortDynamicsMod, only : InitPRTObject + use PRTGenericMod, only : prt_csimpler_allom_hyp use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes @@ -1039,7 +1040,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, c_store) - temp_cohort%laimemory = 0._r8 + temp_cohort%leafmemory = 0._r8 temp_cohort%sapwmemory = 0._r8 temp_cohort%structmemory = 0._r8 cstatus = leaves_on @@ -1048,7 +1049,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & if( prt_params%season_decid(temp_cohort%pft) == itrue .and. & any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf temp_cohort%sapwmemory = c_sapw * stem_drop_fraction temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 @@ -1059,7 +1060,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & if ( prt_params%stress_decid(temp_cohort%pft) == itrue .and. & any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then - temp_cohort%laimemory = c_leaf + temp_cohort%leafmemory = c_leaf temp_cohort%sapwmemory = c_sapw * stem_drop_fraction temp_cohort%structmemory = c_struct * stem_drop_fraction c_leaf = 0._r8 @@ -1140,7 +1141,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & end select select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Equally distribute leaf mass into available age-bins do iage = 1,nleafage @@ -1167,7 +1168,7 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & temp_cohort%coage, temp_cohort%dbh, & - prt_obj, temp_cohort%laimemory,temp_cohort%sapwmemory, temp_cohort%structmemory, & + prt_obj, temp_cohort%leafmemory,temp_cohort%sapwmemory, temp_cohort%structmemory, & cstatus, rstatus, temp_cohort%canopy_trim,temp_cohort%c_area, & 1, csite%spread, bc_in) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 0605767cd6..9d2f3ea9a1 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -100,7 +100,7 @@ module FatesRestartInterfaceMod integer :: ir_coage_co integer :: ir_g_sb_laweight_co integer :: ir_height_co - integer :: ir_laimemory_co + integer :: ir_leafmemory_co integer :: ir_sapwmemory_co integer :: ir_structmemory_co integer :: ir_nplant_co @@ -714,10 +714,10 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - plant height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_height_co ) - call this%set_restart_var(vname='fates_laimemory', vtype=cohort_r8, & + call this%set_restart_var(vname='fates_leafmemory', vtype=cohort_r8, & long_name='ed cohort - target leaf biomass set from prev year', & units='kgC/indiv', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_laimemory_co ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_leafmemory_co ) call this%set_restart_var(vname='fates_sapwmemory', vtype=cohort_r8, & long_name='ed cohort - target sapwood biomass set from prev year', & @@ -1769,7 +1769,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & + rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & @@ -1997,7 +1997,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_dbh_co(io_idx_co) = ccohort%dbh rio_coage_co(io_idx_co) = ccohort%coage rio_height_co(io_idx_co) = ccohort%hite - rio_laimemory_co(io_idx_co) = ccohort%laimemory + rio_leafmemory_co(io_idx_co) = ccohort%leafmemory rio_sapwmemory_co(io_idx_co) = ccohort%sapwmemory rio_structmemory_co(io_idx_co) = ccohort%structmemory rio_g_sb_laweight_co(io_idx_co)= ccohort%g_sb_laweight @@ -2600,7 +2600,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & - rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & + rio_leafmemory_co => this%rvars(ir_leafmemory_co)%r81d, & rio_sapwmemory_co => this%rvars(ir_sapwmemory_co)%r81d, & rio_structmemory_co => this%rvars(ir_structmemory_co)%r81d, & rio_nplant_co => this%rvars(ir_nplant_co)%r81d, & @@ -2802,7 +2802,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%coage = rio_coage_co(io_idx_co) ccohort%g_sb_laweight= rio_g_sb_laweight_co(io_idx_co) ccohort%hite = rio_height_co(io_idx_co) - ccohort%laimemory = rio_laimemory_co(io_idx_co) + ccohort%leafmemory = rio_leafmemory_co(io_idx_co) ccohort%sapwmemory = rio_sapwmemory_co(io_idx_co) ccohort%structmemory= rio_structmemory_co(io_idx_co) ccohort%n = rio_nplant_co(io_idx_co) diff --git a/parteh/PRTAllometricCNPMod.F90 b/parteh/PRTAllometricCNPMod.F90 index 5617d71e5d..3179387590 100644 --- a/parteh/PRTAllometricCNPMod.F90 +++ b/parteh/PRTAllometricCNPMod.F90 @@ -2222,72 +2222,94 @@ end function AllomCNPGrowthDeriv ! ==================================================================================== - subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & - bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, & - grow_leaf,grow_froot,grow_sapw,grow_store) + subroutine TargetAllometryCheck(b0_leaf,b0_fnrt,b0_sapw,b0_store,b0_struct, & + bleaf,bfnrt,bsapw,bstore,bstruct, & + bt_leaf,bt_fnrt,bt_sapw,bt_store,bt_struct, & + carbon_balance,ipft,leaf_status, & + grow_leaf,grow_fnrt,grow_sapw,grow_store,grow_struct) ! Arguments - real(r8),intent(in) :: bleaf !actual - real(r8),intent(in) :: bfroot - real(r8),intent(in) :: bsap + real(r8),intent(in) :: b0_leaf !initial + real(r8),intent(in) :: b0_fnrt + real(r8),intent(in) :: b0_sapw + real(r8),intent(in) :: b0_store + real(r8),intent(in) :: b0_struct + real(r8),intent(in) :: bleaf !actual + real(r8),intent(in) :: bfnrt + real(r8),intent(in) :: bsapw real(r8),intent(in) :: bstore - real(r8),intent(in) :: bdead - real(r8),intent(in) :: bt_leaf !target - real(r8),intent(in) :: bt_froot - real(r8),intent(in) :: bt_sap + real(r8),intent(in) :: bstruct + real(r8),intent(in) :: bt_leaf !target + real(r8),intent(in) :: bt_fnrt + real(r8),intent(in) :: bt_sapw real(r8),intent(in) :: bt_store - real(r8),intent(in) :: bt_dead - logical,intent(out) :: grow_leaf !growth flag - logical,intent(out) :: grow_froot + real(r8),intent(in) :: bt_struct + real(r8),intent(in) :: carbon_balance !remaining carbon balance + integer,intent(in) :: ipft !Plant functional type + integer,intent(in) :: leaf_status !Phenology status + logical,intent(out) :: grow_leaf !growth flag + logical,intent(out) :: grow_fnrt logical,intent(out) :: grow_sapw logical,intent(out) :: grow_store - - if( (bt_leaf - bleaf)>calloc_abs_error) then - write(fates_log(),*) 'leaves are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bleaf,bt_leaf - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bleaf - bt_leaf)>calloc_abs_error) then - ! leaf is above allometry, ignore - grow_leaf = .false. - else - grow_leaf = .true. - end if - - if( (bt_froot - bfroot)>calloc_abs_error) then - write(fates_log(),*) 'fineroots are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bfroot, bt_froot - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bfroot-bt_froot)>calloc_abs_error ) then - grow_froot = .false. - else - grow_froot = .true. - end if - - if( (bt_sap - bsap)>calloc_abs_error) then - write(fates_log(),*) 'sapwood is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bsap, bt_sap - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bsap-bt_sap)>calloc_abs_error ) then - grow_sapw = .false. + logical,intent(out) :: grow_struct + ! Local variables + logical :: fine_leaf + logical :: fine_fnrt + logical :: fine_sapw + logical :: fine_store + logical :: fine_struct + logical :: all_fine + ! Local constants + character(len= 3), parameter :: fmth = '(a)' + character(len=27), parameter :: fmtb = '(a,3(1x,es12.5,1x,a),1x,l1)' + character(len=13), parameter :: fmte = '(a,1x,es12.5)' + character(len=10), parameter :: fmti = '(a,1x,i12)' + + + ! First test whether or not each pool looks reasonable. + fine_leaf = (bt_leaf - bleaf ) <= calloc_abs_error + fine_fnrt = (bt_fnrt - bfnrt ) <= calloc_abs_error + fine_sapw = (bt_sapw - bsapw ) <= calloc_abs_error + fine_store = (bt_store - bstore ) <= calloc_abs_error + fine_struct = (bt_struct - bstruct) <= calloc_abs_error + all_fine = fine_leaf .and. fine_fnrt .and. fine_sapw .and. & + fine_store .and. fine_struct + + ! Decide whether or not to grow tissues (but only if all tissues look fine). + ! We grow only when biomass is less than target biomass (with tolerance). + if (all_fine) then + grow_leaf = ( bleaf - bt_leaf ) <= calloc_abs_error + grow_fnrt = ( bfnrt - bt_fnrt ) <= calloc_abs_error + grow_sapw = ( bsapw - bt_sapw ) <= calloc_abs_error + grow_store = ( bstore - bt_store ) <= calloc_abs_error + grow_struct = ( bstruct - bt_struct ) <= calloc_abs_error else - grow_sapw = .true. - end if - - if( (bt_store - bstore)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bstore,bt_store + ! If anything looks not fine, write a detailed report + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) ' At least one tissue is not on-allometry at the growth step' + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Biomass and on-allometry test (''F'' means problem)' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmth) ' Tissue | Initial | Current | Target | On-allometry' + write(fates_log(),fmt=fmtb) ' Leaf |',b0_leaf ,'|',bleaf ,'|',bt_leaf ,'|',fine_leaf + write(fates_log(),fmt=fmtb) ' Fine root |',b0_fnrt ,'|',bfnrt ,'|',bt_fnrt ,'|',fine_fnrt + write(fates_log(),fmt=fmtb) ' Sap wood |',b0_sapw ,'|',bsapw ,'|',bt_sapw ,'|',fine_sapw + write(fates_log(),fmt=fmtb) ' Storage |',b0_store ,'|',bstore ,'|',bt_store ,'|',fine_store + write(fates_log(),fmt=fmtb) ' Structural |',b0_struct ,'|',bstruct ,'|',bt_struct ,'|',fine_struct + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Ancillary information' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmti) ' PFT = ',ipft + write(fates_log(),fmt=fmti) ' leaf_status = ',leaf_status + write(fates_log(),fmt=fmte) ' carbon_balance = ',carbon_balance + write(fates_log(),fmt=fmte) ' calloc_abs_error = ',calloc_abs_error + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) '======' call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bstore-bt_store)>calloc_abs_error ) then - grow_store = .false. - else - grow_store = .true. end if - if( (bt_dead - bdead)>calloc_abs_error) then - write(fates_log(),*) 'structure not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bdead,bt_dead - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + return end subroutine TargetAllometryCheck diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 5bdf624502..b02c04db94 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -51,6 +51,9 @@ module PRTAllometricCarbonMod use PRTParametersMod , only : prt_params + use EDTypesMod , only : leaves_on + use EDTypesMod , only : leaves_off + implicit none private @@ -114,15 +117,24 @@ module PRTAllometricCarbonMod ! ------------------------------------------------------------------------------------- - type, public, extends(prt_vartypes) :: callom_prt_vartypes + type, public, extends(prt_vartypes) :: callom_prt_vartypes - contains + contains - procedure :: DailyPRT => DailyPRTAllometricCarbon - procedure :: FastPRT => FastPRTAllometricCarbon + procedure :: DailyPRT => DailyPRTAllometricCarbon + procedure :: FastPRT => FastPRTAllometricCarbon end type callom_prt_vartypes - + + type, public, extends(prt_vartypes) :: csimpler_allom_prt_vartypes + + contains + + procedure :: DailyPRT => DailyPRTAllometricCarbonSimpler + procedure :: FastPRT => FastPRTAllometricCarbonSimpler + + end type csimpler_allom_prt_vartypes + ! ------------------------------------------------------------------------------------ ! ! This next class is an extention of the base instance that maps state variables @@ -313,6 +325,9 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: struct_below_target ! dead (structural) biomass below target amount [kgC] real(r8) :: total_below_target ! total biomass below the allometric target [kgC] + real(r8) :: allocation_factor ! allocation factor (relative to demand) to + ! reconstruct tissues + real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC] real(r8) :: store_target_fraction ! ratio between storage and leaf biomass when on allometry [kgC] @@ -338,6 +353,8 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: repro_c0 ! "" real(r8) :: struct_c0 ! "" + logical :: is_deciduous ! Flag to signal this is a deciduous PFT + logical :: grow_struct logical :: grow_leaf ! Are leaves at allometric target and should be grown? logical :: grow_fnrt ! Are fine-roots at allometric target and should be grown? @@ -376,6 +393,10 @@ subroutine DailyPRTAllometricCarbon(this) integer, parameter :: iexp_leaf = 1 ! index 1 is the expanding (i.e. youngest) ! leaf age class, and therefore ! all new allocation goes into that pool + character(len= 9), parameter :: fmti = '(a,1x,i5)' + character(len=13), parameter :: fmt0 = '(a,1x,es12.5)' + character(len=19), parameter :: fmth = '(a,1x,a5,3(1x,a12))' + character(len=22), parameter :: fmtg = '(a,5x,l1,3(1x,es12.5))' real(r8) :: intgr_params(num_bc_in) ! The boundary conditions to this routine, ! are pressed into an array that is also @@ -406,8 +427,11 @@ subroutine DailyPRTAllometricCarbon(this) intgr_params(:) = un_initialized intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval intgr_params(ac_bc_in_id_pft) = real(this%bc_in(ac_bc_in_id_pft)%ival) - - + + + ! Set some logical flags to simplify "if" blocks + is_deciduous = ( prt_params%stress_decid(ipft) == itrue ) .or. & + ( prt_params%season_decid(ipft) == itrue ) nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class @@ -451,11 +475,12 @@ subroutine DailyPRTAllometricCarbon(this) call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) ! Target leaf biomass according to allometry and trimming - if(leaf_status==2) then - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - else - target_leaf_c = 0._r8 - end if + select case (leaf_status) + case (leaves_on) + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + case (leaves_off) + target_leaf_c = 0.0_r8 + end select ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) @@ -464,45 +489,47 @@ subroutine DailyPRTAllometricCarbon(this) call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + ! ----------------------------------------------------------------------------------- ! III. Prioritize some amount of carbon to replace leaf/root turnover - ! Make sure it isnt a negative payment, and either pay what is available + ! Make sure it isn't a negative payment, and either pay what is available ! or forcefully pay from storage. ! ----------------------------------------------------------------------------------- - - if( prt_params%evergreen(ipft) ==1 ) then + if ( is_deciduous ) then + ! Deciduous PFT. Set leaf demand to zero, but keep fine-root demand. + leaf_c_demand = 0.0_r8 + else + ! Evergreen PFT. Try to meet demands for both leaves and fine roots. leaf_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*sum(this%variables(leaf_c_id)%turnover(:))) - else - leaf_c_demand = 0.0_r8 end if - - fnrt_c_demand = max(0.0_r8, & + + fnrt_c_demand = max(0.0_r8, & prt_params%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd)) total_c_demand = leaf_c_demand + fnrt_c_demand - + if (total_c_demand> nearzero ) then ! We pay this even if we don't have the carbon ! Just don't pay so much carbon that storage+carbon_balance can't pay for it + allocation_factor = max(0.0_r8,min(1.0_r8,(store_c+carbon_balance)/total_c_demand)) - leaf_c_flux = min(leaf_c_demand, & - max(0.0_r8,(store_c+carbon_balance)* & - (leaf_c_demand/total_c_demand))) - - ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - ! If we are testing b4b, then we pay this even if we don't have the carbon - fnrt_c_flux = min(fnrt_c_demand, & - max(0.0_r8, (store_c+carbon_balance)* & - (fnrt_c_demand/total_c_demand))) + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_c_demand * allocation_factor + fnrt_c_flux = fnrt_c_demand * allocation_factor - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1) and fine roots + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + ! Remove fluxes from carbon balance. In case we may have drawn carbon from storage, + ! carbon_balance will become negative, in which case we will deplete carbon from + ! storage in the next step. + carbon_balance = carbon_balance - leaf_c_flux - fnrt_c_flux end if ! ----------------------------------------------------------------------------------- @@ -512,18 +539,22 @@ subroutine DailyPRTAllometricCarbon(this) if( carbon_balance < 0.0_r8 ) then + ! Store_c_flux will be negative, so store_c will be depleted store_c_flux = carbon_balance carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux else - - store_below_target = max(target_store_c - store_c,0.0_r8) + ! Accumulate some carbon in storage. If storage is completely depleted, aim to + ! increase storage, but not to replenish completely so we can still use some + ! carbon for growth. + store_below_target = max(0.0_r8,target_store_c - store_c) store_target_fraction = max(0.0_r8, store_c/target_store_c ) store_c_flux = min(store_below_target,carbon_balance * & max(exp(-1.*store_target_fraction**4._r8) - exp( -1.0_r8 ),0.0_r8)) + ! Move carbon from carbon balance to storage carbon_balance = carbon_balance - store_c_flux store_c = store_c + store_c_flux @@ -533,24 +564,28 @@ subroutine DailyPRTAllometricCarbon(this) ! V. If carbon is still available, prioritize some allocation to replace ! the rest of the leaf/fineroot deficit ! carbon balance is guaranteed to be >=0 beyond this point + ! MLO. Renamed demand with below target to make it consistent with the + ! definitions at the variable declaration part. ! ----------------------------------------------------------------------------------- - - leaf_c_demand = max(0.0_r8,(target_leaf_c - sum(leaf_c(1:nleafage)))) - fnrt_c_demand = max(0.0_r8,(target_fnrt_c - fnrt_c)) + leaf_below_target = max(0.0_r8,target_leaf_c - sum(leaf_c(1:nleafage))) + fnrt_below_target = max(0.0_r8,target_fnrt_c - fnrt_c) - total_c_demand = leaf_c_demand + fnrt_c_demand - - if( (carbon_balance > nearzero ) .and. (total_c_demand>nearzero)) then + total_below_target = leaf_below_target + fnrt_below_target + + if ( (carbon_balance > nearzero) .and. (total_below_target > nearzero) ) then + ! Find fraction of carbon to be allocated to leaves and fine roots + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! MLO. Edited the code to switch the order of operations. The previous code would + ! subtract leaf flux from carbon balance before estimating the fine root flux, + ! potentially allowing less fluxes to fine roots than possible. + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor - leaf_c_flux = min(leaf_c_demand, & - carbon_balance*(leaf_c_demand/total_c_demand)) - carbon_balance = carbon_balance - leaf_c_flux leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - fnrt_c_flux = min(fnrt_c_demand, & - carbon_balance*(fnrt_c_demand/total_c_demand)) - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + + carbon_balance = carbon_balance - leaf_c_flux - fnrt_c_flux end if @@ -571,31 +606,23 @@ subroutine DailyPRTAllometricCarbon(this) sapw_below_target + store_below_target if ( total_below_target > nearzero ) then - - if( total_below_target > carbon_balance) then - leaf_c_flux = carbon_balance * leaf_below_target/total_below_target - fnrt_c_flux = carbon_balance * fnrt_below_target/total_below_target - sapw_c_flux = carbon_balance * sapw_below_target/total_below_target - store_c_flux = carbon_balance * store_below_target/total_below_target - else - leaf_c_flux = leaf_below_target - fnrt_c_flux = fnrt_below_target - sapw_c_flux = sapw_below_target - store_c_flux = store_below_target - end if - - carbon_balance = carbon_balance - leaf_c_flux - leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux - - carbon_balance = carbon_balance - fnrt_c_flux - fnrt_c = fnrt_c + fnrt_c_flux - - carbon_balance = carbon_balance - sapw_c_flux - sapw_c = sapw_c + sapw_c_flux - - carbon_balance = carbon_balance - store_c_flux - store_c = store_c + store_c_flux - + ! Find allocation factor based on available carbon and total demand to meet target. + allocation_factor = min(1.0_r8, carbon_balance / total_below_target) + + ! Find fluxes to individual pools + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor + sapw_c_flux = sapw_below_target * allocation_factor + store_c_flux = store_below_target * allocation_factor + + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + sapw_c = sapw_c + sapw_c_flux + store_c = store_c + store_c_flux + + carbon_balance = carbon_balance - & + ( leaf_c_flux + fnrt_c_flux + & + sapw_c_flux + store_c_flux ) end if end if @@ -639,24 +666,15 @@ subroutine DailyPRTAllometricCarbon(this) ! allow actual pools to be above the target, and in these cases, it sends ! a false on the "grow_<>" flag, allowing the plant to grow into these pools. ! It also checks to make sure that structural biomass is not above the target. + ! ( MLO. Removed the check for storage because the same test is done inside + ! sub-routine TargetAllometryCheck.) - if( (target_store_c - store_c)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting' - write(fates_log(),*) 'cbal: ',carbon_balance - write(fates_log(),*) 'near-zero',nearzero - write(fates_log(),*) 'store_c: ',store_c - write(fates_log(),*) 'target c: ',target_store_c - write(fates_log(),*) 'store_c0:', store_c0 - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - - call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, & - store_c, struct_c, & - target_leaf_c, target_fnrt_c, & - target_sapw_c, target_store_c, target_struct_c, & - grow_struct, grow_leaf, grow_fnrt, grow_sapw, grow_store) + call TargetAllometryCheck(sum(leaf_c0(1:nleafage)),fnrt_c0,sapw_c0,store_c0,struct_c0, & + sum(leaf_c(1:nleafage)), fnrt_c, sapw_c,store_c, struct_c, & + target_leaf_c, target_fnrt_c, target_sapw_c, & + target_store_c, target_struct_c, & + carbon_balance,ipft,leaf_status, & + grow_leaf, grow_fnrt, grow_sapw, grow_store, grow_struct) ! -------------------------------------------------------------------------------- ! The numerical integration of growth requires that the instantaneous state @@ -689,19 +707,21 @@ subroutine DailyPRTAllometricCarbon(this) c_pool(repro_c_id) = repro_c c_pool(dbh_id) = dbh - ! Only grow leaves if we are in a "leaf-on" status - if(leaf_status==2) then - c_mask(leaf_c_id) = grow_leaf - else - c_mask(leaf_c_id) = .false. - end if + ! Only grow leaves if we are in a "leaf-on" status. + select case (leaf_status) + case (leaves_on) + c_mask(leaf_c_id) = grow_leaf + case default + c_mask(leaf_c_id) = .false. + end select c_mask(fnrt_c_id) = grow_fnrt c_mask(sapw_c_id) = grow_sapw - c_mask(store_c_id) = grow_store c_mask(struct_c_id) = grow_struct + + c_mask(store_c_id) = grow_store c_mask(repro_c_id) = .true. ! Always calculate reproduction on growth c_mask(dbh_id) = .true. ! Always increment dbh on growth step - + ! When using the Euler method, we keep things simple. We always try ! to make the first integration step to span the entirety of the integration @@ -714,11 +734,12 @@ subroutine DailyPRTAllometricCarbon(this) do_solve_check: do while( ierr .ne. 0 ) deltaC = min(totalC,this%ode_opt_step) - if(ODESolve == 1) then + select_ODESolve: select case (ODESolve) + case (1) call RKF45(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC, & max_trunc_error,intgr_params,c_pool_out,this%ode_opt_step,step_pass) - elseif(ODESolve == 2) then + case (2) call Euler(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC,intgr_params,c_pool_out) ! step_pass = .true. @@ -741,11 +762,11 @@ subroutine DailyPRTAllometricCarbon(this) else this%ode_opt_step = 0.5*deltaC end if - else + case default write(fates_log(),*) 'An integrator was chosen that does not exist' write(fates_log(),*) 'ODESolve = ',ODESolve call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + end select select_ODESolve nsteps = nsteps + 1 @@ -755,17 +776,22 @@ subroutine DailyPRTAllometricCarbon(this) end if if(nsteps > max_substeps ) then - write(fates_log(),*) 'Plant Growth Integrator could not find' - write(fates_log(),*) 'a solution in less than ',max_substeps,' tries' - write(fates_log(),*) 'Aborting' - write(fates_log(),*) 'carbon_balance',carbon_balance - write(fates_log(),*) 'deltaC',deltaC - write(fates_log(),*) 'totalC',totalC - write(fates_log(),*) 'leaf:',grow_leaf,target_leaf_c,target_leaf_c - sum(leaf_c(:)) - write(fates_log(),*) 'fnrt:',grow_fnrt,target_fnrt_c,target_fnrt_c - fnrt_c - write(fates_log(),*) 'sap:',grow_sapw,target_sapw_c, target_sapw_c - sapw_c - write(fates_log(),*) 'store:',grow_store,target_store_c,target_store_c - store_c - write(fates_log(),*) 'dead:',target_struct_c,target_struct_c - struct_c + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=*) 'Plant Growth Integrator could not find' + write(fates_log(),fmt=*) 'a solution in less than ',max_substeps,' tries.' + write(fates_log(),fmt=*) 'Aborting!' + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=fmti) 'Leaf status =',leaf_status + write(fates_log(),fmt=fmt0) 'Carbon_balance =',carbon_balance + write(fates_log(),fmt=fmt0) 'deltaC =',deltaC + write(fates_log(),fmt=fmt0) 'totalC =',totalC + write(fates_log(),fmt=fmth) ' Tissue |', ' Grow',' Current',' Target' ,' Deficit' + write(fates_log(),fmt=fmtg) ' Leaf |', grow_leaf , sum(leaf_c(:)),target_leaf_c , target_leaf_c - sum(leaf_c(:)) + write(fates_log(),fmt=fmtg) ' Fine root |', grow_fnrt , fnrt_c,target_fnrt_c , target_fnrt_c - fnrt_c + write(fates_log(),fmt=fmtg) ' Sapwood |', grow_sapw , sapw_c,target_sapw_c , target_sapw_c - sapw_c + write(fates_log(),fmt=fmtg) ' Storage |', grow_store , store_c,target_store_c , target_store_c - store_c + write(fates_log(),fmt=fmtg) ' Structural |', grow_struct , struct_c,target_struct_c, target_struct_c - struct_c + write(fates_log(),fmt=*) '---~---' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -862,6 +888,562 @@ subroutine DailyPRTAllometricCarbon(this) end subroutine DailyPRTAllometricCarbon ! ===================================================================================== + + + ! ===================================================================================== + + + + + subroutine DailyPRTAllometricCarbonSimpler(this) + + ! ----------------------------------------------------------------------------------- + ! + ! This is the main routine that handles allocation associated with the 1st + ! hypothesis; carbon only, and growth governed by allometry. + ! MLO. This is a slightly simpler alternative I am testing. + ! + ! This routine is explained in the technical documentation in detail. + ! + ! Some points: + ! 1) dbh, while not a PARTEH "state variable", is passed in from FATES (or other + ! model), is integrated along with the mass based state variables, and then + ! passed back to the ecosystem model. It is a "inout" style boundary condition. + ! + ! 2) It is assumed that both growth respiration, and maintenance respiration + ! costs have already been paid, and therefore the "carbon_balance" boundary + ! condition is the net carbon gained by the plant over the coarse of the day. + ! Think of "daily integrated NPP". + ! + ! 3) This routine will completely spend carbon_balance if it enters as a positive + ! value, or replace carbon balance (using storage) if it enters as a negative value. + ! + ! 4) It is assumed that the ecosystem model calling this routine has ensured that + ! the net amount of negative carbon is no greater than that which can be replaced + ! by storage. This routine will crash gracefully if that is not true. + ! + ! 5) Unlike the original sub-routine, here we do not distinguish carbon lost through + ! maintenance from long-term carbon "debt" (i.e. biomass below allometry). We simply + ! seek to bring the plant back to allometry in case there is carbon to do so. We only + ! maintain the priority of replenishing living tissues (plus storage) over structural + ! tissues (which should be sapwood turnover in any case). + ! + ! 6) If there is any carbon balance left after bringing living tissues to allometry, + ! we try to bring heartwood back on allometry. + ! + ! 7) Finally, if carbon is yet still available, it will grow all pools out concurrently + ! including some to reproduction. + ! + ! ---------------------------------------------------------------------------------- + + + ! The class is the only argument + class(csimpler_allom_prt_vartypes) :: this ! this class + + ! ----------------------------------------------------------------------------------- + ! These are local copies of the in/out boundary condition structure + ! ----------------------------------------------------------------------------------- + + real(r8),pointer :: dbh ! Diameter at breast height [cm] + ! this local will point to both in and out bc's + real(r8),pointer :: carbon_balance ! Daily carbon balance for this cohort [kgC] + + real(r8) :: canopy_trim ! The canopy trimming function [0-1] + integer :: ipft ! Plant Functional Type index + + + real(r8) :: target_leaf_c ! target leaf carbon [kgC] + real(r8) :: target_fnrt_c ! target fine-root carbon [kgC] + real(r8) :: target_sapw_c ! target sapwood carbon [kgC] + real(r8) :: target_store_c ! target storage carbon [kgC] + real(r8) :: target_agw_c ! target above ground carbon in woody tissues [kgC] + real(r8) :: target_bgw_c ! target below ground carbon in woody tissues [kgC] + real(r8) :: target_struct_c ! target structural carbon [kgC] + + real(r8) :: sapw_area ! dummy var, x-section area of sapwood [m2] + + real(r8) :: leaf_below_target ! fineroot biomass below target amount [kgC] + real(r8) :: fnrt_below_target ! fineroot biomass below target amount [kgC] + real(r8) :: sapw_below_target ! sapwood biomass below target amount [kgC] + real(r8) :: store_below_target ! storage biomass below target amount [kgC] + real(r8) :: struct_below_target ! dead (structural) biomass below target amount [kgC] + real(r8) :: total_below_target ! total biomass below the allometric target [kgC] + + real(r8) :: allocation_factor ! allocation factor (relative to demand) to + ! reconstruct tissues + + real(r8) :: flux_adj ! adjustment made to growth flux term to minimize error [kgC] + + logical :: step_pass ! Did the integration step pass? + + real(r8) :: leaf_c_flux ! Transfer into leaves at various stages [kgC] + real(r8) :: fnrt_c_flux ! Transfer into fine-roots at various stages [kgC] + real(r8) :: sapw_c_flux ! Transfer into sapwood at various stages [kgC] + real(r8) :: store_c_flux ! Transfer into storage at various stages [kgC] + real(r8) :: repro_c_flux ! Transfer into reproduction at the final stage [kgC] + real(r8) :: struct_c_flux ! Transfer into structure at various stages [kgC] + + real(r8),dimension(max_nleafage) :: leaf_c0 + + ! Initial value of carbon used to determine net flux + real(r8) :: fnrt_c0 ! during this routine + real(r8) :: sapw_c0 ! "" + real(r8) :: store_c0 ! "" + real(r8) :: repro_c0 ! "" + real(r8) :: struct_c0 ! "" + + logical :: is_deciduous ! Flag to signal this is a deciduous PFT + + logical :: grow_struct + logical :: grow_leaf ! Are leaves at allometric target and should be grown? + logical :: grow_fnrt ! Are fine-roots at allometric target and should be grown? + logical :: grow_sapw ! Is sapwood at allometric target and should be grown? + logical :: grow_store ! Is storage at allometric target and should be grown? + + ! integrator variables + real(r8) :: deltaC ! trial value for substep + integer :: ierr ! error flag for allometric growth step + integer :: nsteps ! number of sub-steps + integer :: istep ! current substep index + real(r8) :: totalC ! total carbon allocated over alometric growth step + real(r8) :: hite_out ! dummy height variable + + integer :: i_var ! index for iterating state variables + integer :: i_age ! index for iterating leaf ages + integer :: nleafage ! number of leaf age classifications + integer :: leaf_status ! are leaves on or off? + real(r8) :: leaf_age_flux ! carbon mass flux between leaf age classification pools + + + ! Integrator variables c_pool are "mostly" carbon variables, but c_pool also includes + ! dbh... + ! ----------------------------------------------------------------------------------- + + real(r8),dimension(n_integration_vars) :: c_pool ! Vector of carbon pools passed to integrator + real(r8),dimension(n_integration_vars) :: c_pool_out ! Vector of carbon pools passed back from integrator + logical,dimension(n_integration_vars) :: c_mask ! Mask of active pools during integration + + integer , parameter :: max_substeps = 300 ! Maximum allowable iterations + + real(r8), parameter :: max_trunc_error = 1.0_r8 ! Maximum allowable truncation error + + integer, parameter :: ODESolve = 2 ! 1=RKF45, 2=Euler + + integer, parameter :: iexp_leaf = 1 ! index 1 is the expanding (i.e. youngest) + ! leaf age class, and therefore + ! all new allocation goes into that pool + character(len= 9), parameter :: fmti = '(a,1x,i5)' + character(len=13), parameter :: fmt0 = '(a,1x,es12.5)' + character(len=19), parameter :: fmth = '(a,1x,a5,3(1x,a12))' + character(len=22), parameter :: fmtg = '(a,5x,l1,3(1x,es12.5))' + + real(r8) :: intgr_params(num_bc_in) ! The boundary conditions to this routine, + ! are pressed into an array that is also + ! passed to the integrators + + associate( & + leaf_c => this%variables(leaf_c_id)%val, & + fnrt_c => this%variables(fnrt_c_id)%val(icd), & + sapw_c => this%variables(sapw_c_id)%val(icd), & + store_c => this%variables(store_c_id)%val(icd), & + repro_c => this%variables(repro_c_id)%val(icd), & + struct_c => this%variables(struct_c_id)%val(icd)) + + + ! ----------------------------------------------------------------------------------- + ! 0. + ! Copy the boundary conditions into readable local variables. + ! We don't use pointers for bc's that ar "in" only, only "in-out" and "out" + ! ----------------------------------------------------------------------------------- + + dbh => this%bc_inout(ac_bc_inout_id_dbh)%rval + carbon_balance => this%bc_inout(ac_bc_inout_id_netdc)%rval + + canopy_trim = this%bc_in(ac_bc_in_id_ctrim)%rval + ipft = this%bc_in(ac_bc_in_id_pft)%ival + leaf_status = this%bc_in(ac_bc_in_id_lstat)%ival + + intgr_params(:) = un_initialized + intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval + intgr_params(ac_bc_in_id_pft) = real(this%bc_in(ac_bc_in_id_pft)%ival) + intgr_params(ac_bc_in_id_lstat) = real(this%bc_in(ac_bc_in_id_lstat)%ival,r8) + + ! Set some logical flags to simplify "if" blocks + is_deciduous = ( prt_params%stress_decid(ipft) == itrue ) .or. & + ( prt_params%season_decid(ipft) == itrue ) + + + nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class + + ! ----------------------------------------------------------------------------------- + ! Call the routine that advances leaves in age. + ! This will move a portion of the leaf mass in each + ! age bin, to the next bin. This will not handle movement + ! of mass from the oldest bin into the litter pool, that is something else. + ! ----------------------------------------------------------------------------------- + + call this%AgeLeaves(ipft,sec_per_day) + + ! ----------------------------------------------------------------------------------- + ! I. Remember the values for the state variables at the beginning of this + ! routines. We will then use that to determine their net allocation and reactive + ! transport flux "%net_alloc" at the end. + ! ----------------------------------------------------------------------------------- + + leaf_c0(1:nleafage) = leaf_c(1:nleafage) ! Set initial leaf carbon + fnrt_c0 = fnrt_c ! Set initial fine-root carbon + sapw_c0 = sapw_c ! Set initial sapwood carbon + store_c0 = store_c ! Set initial storage carbon + repro_c0 = repro_c ! Set initial reproductive carbon + struct_c0 = struct_c ! Set initial structural carbon + + + ! ----------------------------------------------------------------------------------- + ! II. Calculate target size of the biomass compartment for a given dbh. + ! ----------------------------------------------------------------------------------- + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + ! Target total above ground biomass in woody/fibrous tissues [kgC] + call bagw_allom(dbh,ipft,target_agw_c) + + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) + + ! Target total dead (structrual) biomass [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + + ! Target storage carbon [kgC,kgC/cm] + call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + + + + ! ----------------------------------------------------------------------------------- + ! II 1/2. Update target biomass based on the leaf status. + ! ----------------------------------------------------------------------------------- + select case (leaf_status) + case (leaves_off) + target_leaf_c = 0.0_r8 + end select + + + ! ----------------------------------------------------------------------------------- + ! III. If carbon is available, bring all the pools as close to the allometry + ! as possible. This also includes the storage pool, even though carbon may + ! be drawn from the storage. + ! ----------------------------------------------------------------------------------- + + ! Identify living organs (plus storage) that are under target. Priority is given + ! to the organs (or storage) that are the most depleted, without a pre-determined + ! sequence. + leaf_below_target = max( 0.0_r8, target_leaf_c - sum(leaf_c(1:nleafage))) + fnrt_below_target = max( 0.0_r8, target_fnrt_c - fnrt_c ) + sapw_below_target = max( 0.0_r8, target_sapw_c - sapw_c ) + store_below_target = max( 0.0_r8, target_store_c - store_c ) + struct_below_target = max( 0.0_r8, target_struct_c - struct_c ) + total_below_target = leaf_below_target + fnrt_below_target + sapw_below_target + & + store_below_target + struct_below_target + + replenish_allom_check: if ( total_below_target > nearzero ) then + ! Available carbon for transfer is the sum of stored carbon and the daily + ! carbon balance. + allocation_factor = min(1.0_r8, (store_c + carbon_balance) / total_below_target ) + + ! Scale flux so pools can be replenished simultaneously. + leaf_c_flux = leaf_below_target * allocation_factor + fnrt_c_flux = fnrt_below_target * allocation_factor + sapw_c_flux = sapw_below_target * allocation_factor + store_c_flux = store_below_target * allocation_factor + struct_c_flux = struct_below_target * allocation_factor + + ! Replenish pools + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + sapw_c = sapw_c + sapw_c_flux + store_c = store_c + store_c_flux + struct_c = struct_c + struct_c_flux + + ! Update carbon balance + carbon_balance = carbon_balance - & + ( leaf_c_flux + fnrt_c_flux + sapw_c_flux + & + store_c_flux + struct_c_flux ) + + end if replenish_allom_check + + + ! ----------------------------------------------------------------------------------- + ! IV. If carbon balance is negative, reduce the storage pool. Otherwise, try to + ! fill the storage pool before growing. + ! ----------------------------------------------------------------------------------- + update_storage: if ( carbon_balance < 0.0_r8 ) then + + + ! If carbon balance is negative, store_c will be depleted. Otherwise, if this is + ! a dormant drought-deciduous plant that somehow managed to score a positive + ! carbon balance with leaves off (very unlikely), then store_c will increase and + ! take all carbon, effectively preventing any chance of growth during lean times. + store_c_flux = carbon_balance + store_c = store_c + store_c_flux + + ! After this operation, carbon_balance should be zero. + carbon_balance = carbon_balance - store_c_flux + + else + ! Non-negative carbon balance. Use left over carbon to fill the storage + ! pool and try to bring it to allometry before trying to grow in size. + store_below_target = max(0.0_r8,target_store_c - store_c) + store_c_flux = min( store_below_target,carbon_balance) + store_c = store_c + store_c_flux + + ! Update carbon balance + carbon_balance = carbon_balance - store_c_flux + + end if update_storage + + + ! ----------------------------------------------------------------------------------- + ! V. If carbon is yet still available ... + ! Our pools are now either on allometry or above (from fusion). + ! We we can increment those pools at or below, + ! including structure and reproduction according to their rates + ! Use an adaptive euler integration. If the error is not nominal, + ! the carbon balance sub-step (deltaC) will be halved and tried again + ! + ! Note that we compare against calloc_abs_error here because it is possible + ! that all the carbon was effectively used up, but a miniscule amount + ! remains due to numerical precision (ie -20 or so), so even though + ! the plant has not been brought to be "on allometry", it thinks it has carbon + ! left to allocate, and thus it must be on allometry when its not. + ! ----------------------------------------------------------------------------------- + if_stature_growth: if ( carbon_balance > calloc_abs_error ) then + + ! This routine checks that actual carbon is not below that targets. It does + ! allow actual pools to be above the target, and in these cases, it sends + ! a false on the "grow_<>" flag, allowing the plant to grow into these pools. + ! It also checks to make sure that structural biomass is not above the target. + call TargetAllometryCheck(sum(leaf_c0(1:nleafage)),fnrt_c0,sapw_c0,store_c0,struct_c0, & + sum(leaf_c(1:nleafage)), fnrt_c, sapw_c,store_c, struct_c, & + target_leaf_c, target_fnrt_c, target_sapw_c, & + target_store_c, target_struct_c, & + carbon_balance,ipft,leaf_status, & + grow_leaf, grow_fnrt, grow_sapw, grow_store, grow_struct) + + ! -------------------------------------------------------------------------------- + ! The numerical integration of growth requires that the instantaneous state + ! variables are passed in as an array. We call it "c_pool". + ! + ! Initialize the adaptive integrator arrays and flags + ! -------------------------------------------------------------------------------- + + ierr = 1 + totalC = carbon_balance + nsteps = 0 + + c_pool(:) = 0.0_r8 ! Zero state variable array + c_mask(:) = .false. ! This mask tells the integrator + ! which indices are active. Its possible + ! that due to fusion, or previous numerical + ! truncation errors, that one of these pools + ! may be larger than its target! We check + ! this, and if true, then we flag that + ! pool to be ignored. c_mask(i) = .false. + ! For grasses, since they don't grow very + ! large and thus won't accumulate such large + ! errors, we always mask as true. + + c_pool(leaf_c_id) = sum(leaf_c(1:nleafage)) + c_pool(fnrt_c_id) = fnrt_c + c_pool(sapw_c_id) = sapw_c + c_pool(store_c_id) = store_c + c_pool(struct_c_id) = struct_c + c_pool(repro_c_id) = repro_c + c_pool(dbh_id) = dbh + + ! Only grow leaves if we are in a "leaf-on" status + if(leaf_status == leaves_on) then + c_mask(leaf_c_id) = grow_leaf + else + c_mask(leaf_c_id) = .false. + end if + c_mask(fnrt_c_id) = grow_fnrt + c_mask(sapw_c_id) = grow_sapw + c_mask(store_c_id) = grow_store + c_mask(struct_c_id) = grow_struct + c_mask(repro_c_id) = .true. ! Always calculate reproduction on growth + c_mask(dbh_id) = .true. ! Always increment dbh on growth step + + + ! When using the Euler method, we keep things simple. We always try + ! to make the first integration step to span the entirety of the integration + ! window for the independent variable (available carbon) + + if(ODESolve == 2) then + this%ode_opt_step = totalC + end if + + do_solve_check: do while( ierr .ne. 0 ) + + deltaC = min(totalC,this%ode_opt_step) + select_ODESolve: select case (ODESolve) + case (1) + call RKF45(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC, & + max_trunc_error,intgr_params,c_pool_out,this%ode_opt_step,step_pass) + + case (2) + call Euler(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC,intgr_params,c_pool_out) + ! step_pass = .true. + + ! When integrating along the allometric curve, we have the luxury of perfect + ! hindsite. Ie, after we have made our step, we can see if the amount + ! of each carbon we have matches the target associated with the new dbh. + ! The following call evaluates how close we are to the allometically defined + ! targets. If we are too far (governed by max_trunc_error), then we + ! pass back the pass/fail flag (step_pass) as false. If false, then + ! we halve the step-size, and then retry. If that step was fine, then + ! we remember the current step size as a good next guess. + + call CheckIntegratedAllometries(c_pool_out(dbh_id),ipft,canopy_trim, & + c_pool_out(leaf_c_id), c_pool_out(fnrt_c_id), c_pool_out(sapw_c_id), & + c_pool_out(store_c_id), c_pool_out(struct_c_id), & + c_mask(leaf_c_id), c_mask(fnrt_c_id), c_mask(sapw_c_id), & + c_mask(store_c_id),c_mask(struct_c_id), max_trunc_error, step_pass) + if(step_pass) then + this%ode_opt_step = deltaC + else + this%ode_opt_step = 0.5*deltaC + end if + case default + write(fates_log(),*) 'An integrator was chosen that does not exist' + write(fates_log(),*) 'ODESolve = ',ODESolve + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select select_ODESolve + + nsteps = nsteps + 1 + + if (step_pass) then ! If true, then step is accepted + totalC = totalC - deltaC + c_pool(:) = c_pool_out(:) + end if + + if(nsteps > max_substeps ) then + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=*) 'Plant Growth Integrator could not find' + write(fates_log(),fmt=*) 'a solution in less than ',max_substeps,' tries.' + write(fates_log(),fmt=*) 'Aborting!' + write(fates_log(),fmt=*) '---~---' + write(fates_log(),fmt=fmti) 'Leaf status =',leaf_status + write(fates_log(),fmt=fmt0) 'Carbon_balance =',carbon_balance + write(fates_log(),fmt=fmt0) 'deltaC =',deltaC + write(fates_log(),fmt=fmt0) 'totalC =',totalC + write(fates_log(),fmt=fmth) ' Tissue |', ' Grow',' Current',' Target' ,' Deficit' + write(fates_log(),fmt=fmtg) ' Leaf |', grow_leaf , sum(leaf_c(:)),target_leaf_c , target_leaf_c - sum(leaf_c(:)) + write(fates_log(),fmt=fmtg) ' Fine root |', grow_fnrt , fnrt_c,target_fnrt_c , target_fnrt_c - fnrt_c + write(fates_log(),fmt=fmtg) ' Sapwood |', grow_sapw , sapw_c,target_sapw_c , target_sapw_c - sapw_c + write(fates_log(),fmt=fmtg) ' Storage |', grow_store , store_c,target_store_c , target_store_c - store_c + write(fates_log(),fmt=fmtg) ' Structural |', grow_struct , struct_c,target_struct_c, target_struct_c - struct_c + write(fates_log(),fmt=*) '---~---' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! + ! TotalC should eventually be whittled down to near zero. + ! The solvers are not perfect, so we can't expect it to be perfectly zero. + ! Note that calloc_abs_error is 1e-9, which is really small (1 microgram of carbon) + ! yet also six orders of magnitude greater than typical rounding errors (~1e-15). + + ! At that point, update the actual states + ! -------------------------------------------------------------------------------- + if_step_pass: if( (totalC < calloc_abs_error) .and. (step_pass) )then + + ierr = 0 + leaf_c_flux = c_pool(leaf_c_id) - sum(leaf_c(1:nleafage)) + fnrt_c_flux = c_pool(fnrt_c_id) - fnrt_c + sapw_c_flux = c_pool(sapw_c_id) - sapw_c + store_c_flux = c_pool(store_c_id) - store_c + struct_c_flux = c_pool(struct_c_id) - struct_c + repro_c_flux = c_pool(repro_c_id) - repro_c + + ! Make an adjustment to flux partitions to make it match remaining c balance + flux_adj = carbon_balance/(leaf_c_flux+fnrt_c_flux+sapw_c_flux + & + store_c_flux+struct_c_flux+repro_c_flux) + + + leaf_c_flux = leaf_c_flux*flux_adj + fnrt_c_flux = fnrt_c_flux*flux_adj + sapw_c_flux = sapw_c_flux*flux_adj + store_c_flux = store_c_flux*flux_adj + struct_c_flux = struct_c_flux*flux_adj + repro_c_flux = repro_c_flux*flux_adj + + carbon_balance = carbon_balance - leaf_c_flux + leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux + + carbon_balance = carbon_balance - fnrt_c_flux + fnrt_c = fnrt_c + fnrt_c_flux + + carbon_balance = carbon_balance - sapw_c_flux + sapw_c = sapw_c + sapw_c_flux + + carbon_balance = carbon_balance - store_c_flux + store_c = store_c + store_c_flux + + carbon_balance = carbon_balance - struct_c_flux + struct_c = struct_c + struct_c_flux + + carbon_balance = carbon_balance - repro_c_flux + repro_c = repro_c + repro_c_flux + + dbh = c_pool(dbh_id) + + if( abs(carbon_balance)>calloc_abs_error ) then + write(fates_log(),*) 'carbon conservation error while integrating pools' + write(fates_log(),*) 'along alometric curve' + write(fates_log(),*) 'carbon_balance = ',carbon_balance,totalC + write(fates_log(),*) 'exiting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end if if_step_pass + + end do do_solve_check + + end if if_stature_growth + + ! Track the net allocations and transport from this routine + ! (the AgeLeaves() routine handled tracking allocation through aging) + + this%variables(leaf_c_id)%net_alloc(icd) = & + this%variables(leaf_c_id)%net_alloc(icd) + (leaf_c(icd) - leaf_c0(icd)) + + this%variables(fnrt_c_id)%net_alloc(icd) = & + this%variables(fnrt_c_id)%net_alloc(icd) + (fnrt_c - fnrt_c0) + + this%variables(sapw_c_id)%net_alloc(icd) = & + this%variables(sapw_c_id)%net_alloc(icd) + (sapw_c - sapw_c0) + + this%variables(store_c_id)%net_alloc(icd) = & + this%variables(store_c_id)%net_alloc(icd) + (store_c - store_c0) + + this%variables(repro_c_id)%net_alloc(icd) = & + this%variables(repro_c_id)%net_alloc(icd) + (repro_c - repro_c0) + + this%variables(struct_c_id)%net_alloc(icd) = & + this%variables(struct_c_id)%net_alloc(icd) + (struct_c - struct_c0) + + + + end associate + + return + end subroutine DailyPRTAllometricCarbonSimpler + + ! ===================================================================================== function AllomCGrowthDeriv(c_pools,c_mask,cbalance,intgr_params) result(dCdx) @@ -1008,78 +1590,92 @@ end function AllomCGrowthDeriv ! ==================================================================================== - subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & - bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, & - grow_dead,grow_leaf,grow_froot,grow_sapw,grow_store) + subroutine TargetAllometryCheck(b0_leaf,b0_fnrt,b0_sapw,b0_store,b0_struct, & + bleaf,bfnrt,bsapw,bstore,bstruct, & + bt_leaf,bt_fnrt,bt_sapw,bt_store,bt_struct, & + carbon_balance,ipft,leaf_status, & + grow_leaf,grow_fnrt,grow_sapw,grow_store,grow_struct) ! Arguments - real(r8),intent(in) :: bleaf !actual - real(r8),intent(in) :: bfroot - real(r8),intent(in) :: bsap + real(r8),intent(in) :: b0_leaf !initial + real(r8),intent(in) :: b0_fnrt + real(r8),intent(in) :: b0_sapw + real(r8),intent(in) :: b0_store + real(r8),intent(in) :: b0_struct + real(r8),intent(in) :: bleaf !actual + real(r8),intent(in) :: bfnrt + real(r8),intent(in) :: bsapw real(r8),intent(in) :: bstore - real(r8),intent(in) :: bdead - real(r8),intent(in) :: bt_leaf !target - real(r8),intent(in) :: bt_froot - real(r8),intent(in) :: bt_sap + real(r8),intent(in) :: bstruct + real(r8),intent(in) :: bt_leaf !target + real(r8),intent(in) :: bt_fnrt + real(r8),intent(in) :: bt_sapw real(r8),intent(in) :: bt_store - real(r8),intent(in) :: bt_dead - logical,intent(out) :: grow_leaf !growth flag - logical,intent(out) :: grow_froot + real(r8),intent(in) :: bt_struct + real(r8),intent(in) :: carbon_balance !remaining carbon balance + integer,intent(in) :: ipft !Plant functional type + integer,intent(in) :: leaf_status !Phenology status + logical,intent(out) :: grow_leaf !growth flag + logical,intent(out) :: grow_fnrt logical,intent(out) :: grow_sapw logical,intent(out) :: grow_store - logical,intent(out) :: grow_dead - - if( (bt_leaf - bleaf)>calloc_abs_error) then - write(fates_log(),*) 'leaves are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bleaf,bt_leaf - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bleaf - bt_leaf)>calloc_abs_error) then - ! leaf is above allometry, ignore - grow_leaf = .false. - else - grow_leaf = .true. - end if - - if( (bt_froot - bfroot)>calloc_abs_error) then - write(fates_log(),*) 'fineroots are not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bfroot, bt_froot - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bfroot-bt_froot)>calloc_abs_error ) then - grow_froot = .false. - else - grow_froot = .true. - end if - - if( (bt_sap - bsap)>calloc_abs_error) then - write(fates_log(),*) 'sapwood is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bsap, bt_sap - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bsap-bt_sap)>calloc_abs_error ) then - grow_sapw = .false. + logical,intent(out) :: grow_struct + ! Local variables + logical :: fine_leaf + logical :: fine_fnrt + logical :: fine_sapw + logical :: fine_store + logical :: fine_struct + logical :: all_fine + ! Local constants + character(len= 3), parameter :: fmth = '(a)' + character(len=27), parameter :: fmtb = '(a,3(1x,es12.5,1x,a),1x,l1)' + character(len=13), parameter :: fmte = '(a,1x,es12.5)' + character(len=10), parameter :: fmti = '(a,1x,i12)' + + + ! First test whether or not each pool looks reasonable. + fine_leaf = (bt_leaf - bleaf ) <= calloc_abs_error + fine_fnrt = (bt_fnrt - bfnrt ) <= calloc_abs_error + fine_sapw = (bt_sapw - bsapw ) <= calloc_abs_error + fine_store = (bt_store - bstore ) <= calloc_abs_error + fine_struct = (bt_struct - bstruct) <= calloc_abs_error + all_fine = fine_leaf .and. fine_fnrt .and. fine_sapw .and. & + fine_store .and. fine_struct + + ! Decide whether or not to grow tissues (but only if all tissues look fine). + ! We grow only when biomass is less than target biomass (with tolerance). + if (all_fine) then + grow_leaf = ( bleaf - bt_leaf ) <= calloc_abs_error + grow_fnrt = ( bfnrt - bt_fnrt ) <= calloc_abs_error + grow_sapw = ( bsapw - bt_sapw ) <= calloc_abs_error + grow_store = ( bstore - bt_store ) <= calloc_abs_error + grow_struct = ( bstruct - bt_struct ) <= calloc_abs_error else - grow_sapw = .true. - end if - - if( (bt_store - bstore)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bstore,bt_store + ! If anything looks not fine, write a detailed report + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) ' At least one tissue is not on-allometry at the growth step' + write(fates_log(),fmt=fmth) '======' + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Biomass and on-allometry test (''F'' means problem)' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmth) ' Tissue | Initial | Current | Target | On-allometry' + write(fates_log(),fmt=fmtb) ' Leaf |',b0_leaf ,'|',bleaf ,'|',bt_leaf ,'|',fine_leaf + write(fates_log(),fmt=fmtb) ' Fine root |',b0_fnrt ,'|',bfnrt ,'|',bt_fnrt ,'|',fine_fnrt + write(fates_log(),fmt=fmtb) ' Sap wood |',b0_sapw ,'|',bsapw ,'|',bt_sapw ,'|',fine_sapw + write(fates_log(),fmt=fmtb) ' Storage |',b0_store ,'|',bstore ,'|',bt_store ,'|',fine_store + write(fates_log(),fmt=fmtb) ' Structural |',b0_struct ,'|',bstruct ,'|',bt_struct ,'|',fine_struct + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) ' Ancillary information' + write(fates_log(),fmt=fmth) '------' + write(fates_log(),fmt=fmti) ' PFT = ',ipft + write(fates_log(),fmt=fmti) ' leaf_status = ',leaf_status + write(fates_log(),fmt=fmte) ' carbon_balance = ',carbon_balance + write(fates_log(),fmt=fmte) ' calloc_abs_error = ',calloc_abs_error + write(fates_log(),fmt=fmth) '' + write(fates_log(),fmt=fmth) '======' call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( ( bstore-bt_store)>calloc_abs_error ) then - grow_store = .false. - else - grow_store = .true. end if - - if( (bt_dead - bdead)>calloc_abs_error) then - write(fates_log(),*) 'structure not on-allometry at the growth step' - write(fates_log(),*) 'exiting',bdead,bt_dead - call endrun(msg=errMsg(sourcefile, __LINE__)) - elseif( (bdead-bt_dead)> calloc_abs_error) then - grow_dead = .false. - else - grow_dead = .true. - end if - return end subroutine TargetAllometryCheck @@ -1100,5 +1696,20 @@ subroutine FastPRTAllometricCarbon(this) end subroutine FastPRTAllometricCarbon + ! ===================================================================================== + + subroutine FastPRTAllometricCarbonSimpler(this) + + implicit none + class(csimpler_allom_prt_vartypes) :: this ! this class + + ! This routine does nothing, because in the carbon only allometric RT model + ! we currently don't have any fast-timestep processes + ! Think of this as a stub. + + + return + end subroutine FastPRTAllometricCarbonSimpler + end module PRTAllometricCarbonMod diff --git a/parteh/PRTGenericMod.F90 b/parteh/PRTGenericMod.F90 index 3dab9563a3..4d686d8ac6 100644 --- a/parteh/PRTGenericMod.F90 +++ b/parteh/PRTGenericMod.F90 @@ -66,6 +66,7 @@ module PRTGenericMod ! These should each have their own module ! ------------------------------------------------------------------------------------- + integer, parameter, public :: prt_csimpler_allom_hyp = 0 integer, parameter, public :: prt_carbon_allom_hyp = 1 integer, parameter, public :: prt_cnp_flex_allom_hyp = 2 ! Still under development diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90 index 04a0f5dda0..d6884ab468 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -100,7 +100,8 @@ module PRTParametersMod real(r8), allocatable :: c2b(:) ! Carbon to biomass multiplier [kg/kgC] real(r8), allocatable :: wood_density(:) ! wood density g cm^-3 ... - real(r8), allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) + + integer , allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) real(r8), allocatable :: crown(:) ! fraction of the height of the plant ! that is occupied by crown real(r8), allocatable :: slamax(:) ! Maximum specific leaf area of plant (at bottom) [m2/gC] diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 208ff848fb..cbc4df5d28 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -17,7 +17,9 @@ module PRTInitParamsFatesMod use FatesGlobals, only : fates_log use shr_log_mod, only : errMsg => shr_log_errMsg use EDPftvarcon, only : EDPftvarcon_inst - use PRTGenericMod, only : prt_cnp_flex_allom_hyp,prt_carbon_allom_hyp + use PRTGenericMod, only : prt_cnp_flex_allom_hyp + use PRTGenericMod, only : prt_carbon_allom_hyp + use PRTGenericMod, only : prt_csimpler_allom_hyp use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : bagw_allom @@ -441,8 +443,11 @@ subroutine PRTReceivePFT(fates_params) name = 'fates_woody' call fates_params%RetreiveParameterAllocate(name=name, & - data=prt_params%woody) - + data=tmpreal) + allocate(prt_params%woody(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%woody) + deallocate(tmpreal) + name = 'fates_wood_density' call fates_params%RetreiveParameterAllocate(name=name, & data=prt_params%wood_density) @@ -1014,12 +1019,12 @@ subroutine PRTCheckParams(is_master) ! Check to make sure the organ ids are valid if this is the ! cnp_flex_allom_hypothesis - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do io = 1,norgans if(prt_params%organ_id(io) == repro_organ) then - write(fates_log(),*) 'with flexible cnp or c-only alloc hypothesese' + write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'reproductive tissues are a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' @@ -1036,7 +1041,7 @@ subroutine PRTCheckParams(is_master) end if end do - end if + end select pftloop: do ipft = 1,npft @@ -1076,13 +1081,13 @@ subroutine PRTCheckParams(is_master) ! Check if woody plants have a structural biomass (agb) intercept ! ---------------------------------------------------------------------------------- if ( ( prt_params%allom_agb1(ipft) <= tiny(prt_params%allom_agb1(ipft)) ) .and. & - ( int(prt_params%woody(ipft)) .eq. 1 ) ) then + ( prt_params%woody(ipft) .eq. 1 ) ) then write(fates_log(),*) 'Woody plants are expected to have a non-zero intercept' write(fates_log(),*) ' in the diameter to AGB allometry equations' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' allom_agb1: ',prt_params%allom_agb1(ipft) - write(fates_log(),*) ' woody: ',int(prt_params%woody(ipft)) + write(fates_log(),*) ' woody: ',prt_params%woody(ipft) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -1091,7 +1096,7 @@ subroutine PRTCheckParams(is_master) ! Check if non-woody plants have structural biomass (agb) intercept ! ---------------------------------------------------------------------------------- ! if ( ( prt_params%allom_agb1(ipft) > tiny(prt_params%allom_agb1(ipft)) ) .and. & -! ( int(prt_params%woody(ipft)) .ne. 1 ) ) then +! ( prt_params%woody(ipft) .ne. 1 ) ) then ! ! write(fates_log(),*) 'Non-woody plants are expected to have a zero intercept' ! write(fates_log(),*) ' in the diameter to AGB allometry equations' @@ -1100,7 +1105,7 @@ subroutine PRTCheckParams(is_master) ! write(fates_log(),*) ' woody tissues (sap and structural dead wood).' ! write(fates_log(),*) ' PFT#: ',ipft ! write(fates_log(),*) ' allom_agb1: ',prt_params%allom_agb1(ipft) -! write(fates_log(),*) ' woody: ',int(prt_params%woody(ipft)) +! write(fates_log(),*) ' woody: ',prt_params%woody(ipft) ! write(fates_log(),*) ' Aborting' ! call endrun(msg=errMsg(sourcefile, __LINE__)) ! @@ -1128,9 +1133,8 @@ subroutine PRTCheckParams(is_master) ! should not be re-translocating mass upon turnover. ! Note to advanced users. Feel free to remove these checks... ! ------------------------------------------------------------------- - - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) do i = 1,norgans io = prt_params%organ_id(i) @@ -1165,10 +1169,11 @@ subroutine PRTCheckParams(is_master) end if end do - end if - - if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then - + end select + + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) + ! Make sure nutrient storage fractions are positive if( prt_params%nitr_store_ratio(ipft) < 0._r8 ) then write(fates_log(),*) 'With parteh allometric CNP hypothesis' @@ -1235,9 +1240,8 @@ subroutine PRTCheckParams(is_master) end if end do - - end if - + + end select ! Growth respiration ! if (parteh_mode .eq. prt_carbon_allom_hyp) then @@ -1258,8 +1262,8 @@ subroutine PRTCheckParams(is_master) ! end if ! end if - if ((hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) .or. & - (hlm_parteh_mode .eq. prt_carbon_allom_hyp) ) then + select case (hlm_parteh_mode) + case (prt_csimpler_allom_hyp,prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp) ! The first nitrogen stoichiometry is used in all cases if ( (any(prt_params%nitr_stoich_p1(ipft,:) < 0.0_r8)) .or. & (any(prt_params%nitr_stoich_p1(ipft,:) >= 1.0_r8))) then @@ -1269,9 +1273,10 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select - if(hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then + select case (hlm_parteh_mode) + case (prt_cnp_flex_allom_hyp) do i = 1,norgans if ( (prt_params%nitr_stoich_p1(ipft,i) < 0._r8) .or. & @@ -1309,7 +1314,7 @@ subroutine PRTCheckParams(is_master) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end select ! Check turnover time-scales