diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 00a969a78c..3b8d94e6ed 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1067,7 +1067,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort%c_area/currentPatch%total_canopy_area) currentPatch%layer_height_profile(L,ft,iv) = currentPatch%layer_height_profile(L,ft,iv) + (remainder * fleaf * & currentCohort%c_area/currentPatch%total_canopy_area*(layer_top_hite+layer_bottom_hite)/2.0_r8) - write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(L,ft,iv) + if ( DEBUG ) write(fates_log(), *) 'LHP', currentPatch%layer_height_profile(L,ft,iv) if(currentCohort%dbh <= 0._r8.or.currentCohort%n == 0._r8)then write(fates_log(), *) 'ED: dbh or n is zero in clmedlink', currentCohort%dbh,currentCohort%n endif diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 199682bdb7..2237553ccd 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -186,6 +186,9 @@ subroutine allocate_live_biomass(cc_p,mode) ! accounts for the fact that live biomass may decline in the off-season, ! making leaf_memory unrealistic. real(r8) :: ratio_balive ! ratio between root+shoot biomass now and root+shoot biomass when leaves fell off. + real(r8) :: new_bl + real(r8) :: new_br + real(r8) :: new_bsw integer :: ft ! functional type integer :: leaves_off_switch @@ -220,68 +223,69 @@ subroutine allocate_live_biomass(cc_p,mode) ! Use different proportions if the leaves are on vs off if(leaves_off_switch==0)then - ! Tracking npp/gpp diagnostics only occur after growth derivatives is called - if(mode==1)then - ! it will not be able to put out as many leaves as it had previous timestep - currentcohort%npp_leaf = currentcohort%npp_leaf + & - max(0.0_r8,currentcohort%balive*leaf_frac - currentcohort%bl)/freq_day - end if - currentcohort%bl = currentcohort%balive*leaf_frac + new_bl = currentcohort%balive*leaf_frac + + new_br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac + + new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & + currentcohort%laimemory)*leaf_frac !diagnose the root and stem biomass from the functional balance hypothesis. This is used when the leaves are !fully on. - if(mode==1)then + if(mode==1)then + currentcohort%npp_leaf = currentcohort%npp_leaf + & + max(0.0_r8,new_bl - currentcohort%bl) / freq_day + currentcohort%npp_froot = currentcohort%npp_froot + & - max(0._r8,pftcon%froot_leaf(ft)*(currentcohort%balive+currentcohort%laimemory)*leaf_frac - currentcohort%br) / & - freq_day + max(0._r8,new_br - currentcohort%br) / freq_day - currentcohort%npp_bsw = max(0._r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & - currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/freq_day + currentcohort%npp_bsw = max(0.0_r8, new_bsw - currentcohort%bsw)/freq_day currentcohort%npp_bdead = currentCohort%dbdeaddt end if + + currentcohort%bl = new_bl + currentcohort%br = new_br + currentcohort%bsw = new_bsw - currentcohort%br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac - currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & - currentcohort%laimemory)*leaf_frac - - else ! Leaves are on (leaves_off_switch==1) - - !the purpose of this section is to figure out the root and stem biomass when the leaves are off - !at this point, we know the former leaf mass (laimemory) and the current alive mass - !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted - !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass - !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF + else ! Leaves are off (leaves_off_switch==1) - currentcohort%bl = 0.0_r8 + !the purpose of this section is to figure out the root and stem biomass when the leaves are off + !at this point, we know the former leaf mass (laimemory) and the current alive mass + !because balive may decline in the off-season, we need to adjust the + !root and stem biomass that are predicted from the laimemory, for the fact that we now might + !not have enough live biomass to support the hypothesized root mass + !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF + ideal_balive = currentcohort%laimemory * pftcon%froot_leaf(ft) + & currentcohort%laimemory* EDecophyscon%sapwood_ratio(ft) * currentcohort%hite - currentcohort%br = pftcon%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * leaf_frac - currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + & - currentcohort%laimemory)*leaf_frac - - ratio_balive = currentcohort%balive / ideal_balive - currentcohort%br = currentcohort%br * ratio_balive - currentcohort%bsw = currentcohort%bsw * ratio_balive + ratio_balive = currentcohort%balive / ideal_balive + + new_br = pftcon%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * & + leaf_frac * ratio_balive + new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite * & + (ideal_balive + currentcohort%laimemory) * leaf_frac * ratio_balive ! Diagnostics if(mode==1)then currentcohort%npp_froot = currentcohort%npp_froot + & - max(0.0_r8,pftcon%froot_leaf(ft)*(ideal_balive + & - currentcohort%laimemory)*leaf_frac*ratio_balive-currentcohort%br)/freq_day + max(0.0_r8,new_br-currentcohort%br)/freq_day - currentcohort%npp_bsw = max(0.0_r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + & - currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/freq_day + currentcohort%npp_bsw = max(0.0_r8, new_bsw-currentcohort%bsw)/freq_day currentcohort%npp_bdead = currentCohort%dbdeaddt end if + currentcohort%bl = 0.0_r8 + currentcohort%br = new_br + currentcohort%bsw = new_bsw + endif if (abs(currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw)>1e-12) then @@ -1041,12 +1045,11 @@ subroutine copy_cohort( currentCohort,copyc ) n%npp_acc_hold = o%npp_acc_hold n%npp_tstep = o%npp_tstep + n%npp_acc = o%npp_acc if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ia ',o%npp_acc if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ib ',o%resp_acc - n%npp_acc = o%npp_acc - n%resp_tstep = o%resp_tstep n%resp_acc = o%resp_acc n%resp_acc_hold = o%resp_acc_hold diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index a496767038..430da23a4b 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -21,6 +21,12 @@ module EDPhysiologyMod use EDTypesMod , only : ncwd, cp_nlevcan, numpft_ed, senes use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type + use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun + use FatesGlobals , only : fates_log + + + implicit none private @@ -40,7 +46,8 @@ module EDPhysiologyMod public :: flux_into_litter_pools logical, parameter :: DEBUG = .false. ! local debug flag - + character(len=*), parameter, private :: sourcefile = & + __FILE__ ! ============================================================================ contains @@ -817,12 +824,12 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) currentCohort%carbon_balance = currentCohort%npp_acc_hold - & currentCohort%md * EDecophyscon%leaf_stor_priority(currentCohort%pft) - ! Allowing only carbon from NPP pool to account for npp flux into the maintenance pools + ! Allowing only carbon from NPP pool to account for npp flux into the maintenance turnover pools ! ie this does not include any use of storage carbon or balive to make up for missing carbon balance in the transfer - currentCohort%npp_leaf = min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, & - currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)) - currentCohort%npp_froot = min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, & - currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft)) + currentCohort%npp_leaf = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%leaf_md/currentCohort%md, & + currentCohort%leaf_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))) + currentCohort%npp_froot = max(0.0_r8,min(currentCohort%npp_acc_hold*currentCohort%root_md/currentCohort%md, & + currentCohort%root_md*EDecophyscon%leaf_stor_priority(currentCohort%pft))) if (Bleaf(currentCohort) > 0._r8)then @@ -839,21 +846,26 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) !what is the flux into the store? currentCohort%storage_flux = currentCohort%carbon_balance * f_store + currentCohort%npp_store = currentCohort%carbon_balance * f_store if ( DEBUG ) write(fates_log(),*) 'EDphys B ',f_store !what is the tax on the carbon available for growth? currentCohort%carbon_balance = currentCohort%carbon_balance * (1.0_r8 - f_store) else !cbalance is negative. Take C out of store to pay for maintenance respn. + currentCohort%storage_flux = currentCohort%carbon_balance + + ! Note that npp_store only tracks the flux between NPP and storage. Storage can + ! also be drawn down to support some turnover demand. + currentCohort%npp_store = min(0.0_r8,currentCohort%npp_acc_hold) + currentCohort%carbon_balance = 0._r8 endif else - currentCohort%storage_flux = 0._r8 - currentCohort%carbon_balance = 0._r8 - write(fates_log(),*) 'ED: no leaf area in gd',currentCohort%n,currentCohort%bdead, & - currentCohort%dbh,currentCohort%balive + write(fates_log(),*) 'No target leaf area in GrowthDerivs? Bleaf(cohort) <= 0?' + call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -947,7 +959,6 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) endif currentCohort%npp_bseed = currentCohort%seed_prod - currentCohort%npp_store = max(0.0_r8,currentCohort%storage_flux) ! calculate change in diameter and height currentCohort%ddbhdt = currentCohort%dbdeaddt * dDbhdBd(currentCohort) diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index 78563a3a2d..bd2437c92f 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -9,14 +9,18 @@ module EDAccumulateFluxesMod ! Rosie Fisher. March 2014. ! ! !USES: + use abortutils, only : endrun + use shr_log_mod , only : errMsg => shr_log_errMsg implicit none private ! public :: AccumulateFluxes_ED logical :: DEBUG = .false. ! for debugging this module - !------------------------------------------------------------------------------ + character(len=*), parameter, private :: sourcefile = & + __FILE__ + contains !------------------------------------------------------------------------------ @@ -33,6 +37,8 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) use EDTypesMod , only : ed_patch_type, ed_cohort_type, & ed_site_type, AREA use FatesInterfaceMod , only : bc_in_type,bc_out_type + use, intrinsic :: IEEE_ARITHMETIC + ! ! !ARGUMENTS integer, intent(in) :: nsites @@ -67,12 +73,11 @@ subroutine AccumulateFluxes_ED(nsites, sites, bc_in, bc_out, dt_time) ! _tstep fluxes are KgC/indiv/timestep _acc are KgC/indiv/day if ( DEBUG ) then - write(iulog,*) 'EDAccumFlux 64 ',ccohort%npp_acc, & - ccohort%npp_tstep + write(iulog,*) 'EDAccumFlux 64 ',ccohort%npp_tstep write(iulog,*) 'EDAccumFlux 66 ',ccohort%gpp_tstep write(iulog,*) 'EDAccumFlux 67 ',ccohort%resp_tstep endif - + ccohort%npp_acc = ccohort%npp_acc + ccohort%npp_tstep ccohort%gpp_acc = ccohort%gpp_acc + ccohort%gpp_tstep ccohort%resp_acc = ccohort%resp_acc + ccohort%resp_tstep diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 8ac4a51b36..8283a4d52d 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -192,7 +192,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) end if enddo enddo - + !weight patch level output BTRAN for the bc_out(s)%btran_pa(ifp) = 0.0_r8 do ft = 1,numpft_ed @@ -203,10 +203,11 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out) bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft_ed end if enddo - - temprootr = sum(bc_out(s)%rootr_pagl(ifp,:)) + + ! While the in-pft root profiles summed to unity, averaging them weighted + ! by conductance, or not, will break sum to unity. Thus, re-normalize. + temprootr = sum(bc_out(s)%rootr_pagl(ifp,1:cp_numlevgrnd)) if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then - write(iulog,*) 'error with rootr in canopy fluxes',temprootr,sum(pftgs),sum(cpatch%rootr_ft(1:2,:),dim=2) do j = 1,cp_numlevgrnd bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr enddo diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index 91d8d844fd..7a9e3e7467 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -206,7 +206,6 @@ subroutine FATES_BGC_Carbon_Balancecheck(nsites, sites, bc_in, is_beg_day, dtime sites(s)%fire_c_to_atm*SHR_CONST_CDAY) sites(s)%cbal_err_fates = sites(s)%cbal_err_fates / SHR_CONST_CDAY - sites(s)%cbal_err_bgc = sites(s)%totbgcc - & sites(s)%totbgcc_old - & (sites(s)%fates_to_bgc_last_ts*SHR_CONST_CDAY - & diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index b7bf5edb9b..bf1f7e562f 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -58,10 +58,16 @@ module FatesConstantsMod real(fates_r8), parameter :: t_water_freeze_k_triple = 273.16_fates_r8 + ! For numerical inquiry + real(fates_r8), parameter :: fates_huge = huge(g_per_kg) + + real(fates_r8), parameter :: fates_tiny = tiny(g_per_kg) + ! Geometric Constants ! PI real(fates_r8), parameter :: pi_const = 3.14159265359_fates_r8 + end module FatesConstantsMod