Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixes to carbon accounting and lr2 machine files #174

Merged
merged 10 commits into from
Feb 6, 2017
Merged
4 changes: 0 additions & 4 deletions cime/cime_config/cesm/machines/config_machines.xml
Original file line number Diff line number Diff line change
Expand Up @@ -861,10 +861,6 @@
</batch_system>
<mpirun mpilib="mpi-serial">
<executable></executable>
<arguments>
<arg name="num_tasks">-np {{ num_tasks }}</arg>
<arg name="tasks_per_node"> -npernode {{ tasks_per_node }}</arg>
</arguments>
</mpirun>
<mpirun mpilib="default">
<executable>mpirun</executable>
Expand Down
2 changes: 1 addition & 1 deletion components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
77 changes: 40 additions & 37 deletions components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
33 changes: 22 additions & 11 deletions components/clm/src/ED/biogeochem/EDPhysiologyMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand Down
13 changes: 9 additions & 4 deletions components/clm/src/ED/biogeophys/EDAccumulateFluxesMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

!------------------------------------------------------------------------------
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
9 changes: 5 additions & 4 deletions components/clm/src/ED/biogeophys/EDBtranMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
1 change: 0 additions & 1 deletion components/clm/src/ED/main/ChecksBalancesMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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 - &
Expand Down
6 changes: 6 additions & 0 deletions components/clm/src/ED/main/FatesConstantsMod.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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