Skip to content

Commit

Permalink
Merge branch 'peterdschwartz/lnd/fix-balance-errors' into next (PR #6235
Browse files Browse the repository at this point in the history 
)

Carbon balance was never checked correctly and let through negative values.
Initialization of Ecosystem variables needed to be re-adjusted, and balance
checking for FATES needed to be changed, which causes round-off differences
in CMASS_BALANCE_ERROR for fate_cold_allvars test.

Also included divide-by-zero check in PhenologyMod.F90 to fix fpe in debug mode.

[non-BFB] for one FATES test.
Fixes #6120
Fixes #6203
Fixes #6177
  • Loading branch information
@bishtgautam
bishtgautam committed May 30, 2024
2 parents 3e356e8 + e6097d2 commit 196a9b2
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 62 deletions.
106 changes: 65 additions & 41 deletions components/elm/src/biogeochem/EcosystemBalanceCheckMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,44 +285,41 @@ subroutine ColCBalanceCheck(bounds, &
end do ! end of columns loop

! Consider adapting this check to be fates compliant (rgk 04-2017)
if (.not. use_fates) then
#ifndef _OPENACC
if (err_found .and. nstep > 1) then
c = err_index
write(iulog,*)'column cbalance error = ', col_errcb(c), c
write(iulog,*)'Latdeg,Londeg = ',grc_pp%latdeg(col_pp%gridcell(c)),grc_pp%londeg(col_pp%gridcell(c))
write(iulog,*)'input = ',col_cinputs(c)*dt
write(iulog,*)'output = ',col_coutputs(c)*dt
write(iulog,*)'er = ',er(c)*dt,col_cf%hr(c)*dt
write(iulog,*)'fire = ',col_fire_closs(c)*dt
write(iulog,*)'hrv_to_atm = ',col_hrv_xsmrpool_to_atm(c)*dt
write(iulog,*)'leach = ',som_c_leached(c)*dt
write(iulog,*)'begcb = ',col_begcb(c)
write(iulog,*)'endcb = ',col_endcb(c),col_cs%totsomc(c)
write(iulog,*)'totsomc = ',col_cs%totsomc(c)
write(iulog,*)'delta store = ',col_endcb(c)-col_begcb(c)

if (ero_ccycle) then
write(iulog,*)'erosion = ',som_c_yield(c)*dt
end if

if (use_pflotran .and. pf_cmode) then
write(iulog,*)'pf_delta_decompc = ',col_decompc_delta(c)*dt
end if

if (use_pflotran .and. pf_cmode) then
write(iulog,*)'pf_delta_decompc = ',col_decompc_delta(c)*dt
end if

call endrun(msg=errMsg(__FILE__, __LINE__))
else
if (masterproc .and. nstep < 2) then
write(iulog,*) '--WARNING-- skipping CN balance check for first timestep'
end if
end if
if (err_found .and. nstep > 1) then
c = err_index
write(iulog,*)'column cbalance error = ', col_errcb(c), c
write(iulog,*)'Latdeg,Londeg = ',grc_pp%latdeg(col_pp%gridcell(c)),grc_pp%londeg(col_pp%gridcell(c))
write(iulog,*)'input = ',col_cinputs(c)*dt
write(iulog,*)'output = ',col_coutputs(c)*dt
write(iulog,*)'er = ',er(c)*dt,col_cf%hr(c)*dt
write(iulog,*)'fire = ',col_fire_closs(c)*dt
write(iulog,*)'hrv_to_atm = ',col_hrv_xsmrpool_to_atm(c)*dt
write(iulog,*)'leach = ',som_c_leached(c)*dt
write(iulog,*)'begcb = ',col_begcb(c)
write(iulog,*)'endcb = ',col_endcb(c),col_cs%totsomc(c)
write(iulog,*)'totsomc = ',col_cs%totsomc(c)
write(iulog,*)'delta store = ',col_endcb(c)-col_begcb(c)

if (ero_ccycle) then
write(iulog,*)'erosion = ',som_c_yield(c)*dt
end if

if (use_pflotran .and. pf_cmode) then
write(iulog,*)'pf_delta_decompc = ',col_decompc_delta(c)*dt
end if

if (use_pflotran .and. pf_cmode) then
write(iulog,*)'pf_delta_decompc = ',col_decompc_delta(c)*dt
end if

call endrun(msg=errMsg(__FILE__, __LINE__))
else
if (masterproc .and. nstep < 2) then
write(iulog,*) '--WARNING-- skipping CN balance check for first timestep'
end if
end if
#endif
end if


end associate

Expand Down Expand Up @@ -624,7 +621,6 @@ subroutine ColPBalanceCheck(bounds, &
! set time steps
dt = dtime_mod
kyr = year_curr; kmo = mon_curr; kda = day_curr; mcsec = secs_curr;

err_found = .false.

if(.not.use_fates)then
Expand Down Expand Up @@ -907,7 +903,14 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
grc_som_c_leached => grc_cf%som_c_leached , & ! Output: [real(r8) (:) ] (gC/m^2/s)total SOM C loss from vertical transport
grc_som_c_yield => grc_cf%somc_yield , & ! Output: [real(r8) (:) ] (gC/m^2/s)total SOM C loss by erosion
grc_cinputs => grc_cf%cinputs , & ! Output: [real(r8) (:) ] (gC/m2/s) column-level C inputs
grc_coutputs => grc_cf%coutputs & ! Output: [real(r8) (:) ] (gC/m2/s) column-level C outputs
grc_coutputs => grc_cf%coutputs , & ! Output: [real(r8) (:) ] (gC/m2/s) column-level C outputs
beg_totpftc => grc_cs%beg_totpftc , & ! Input: [real(r8) (:)] (gC/m2) patch-level carbon aggregated to column level, incl veg and cpool
beg_cwdc => grc_cs%beg_cwdc , & ! Input: [real(r8) (:)] (gC/m2) total column coarse woody debris carbon
beg_totsomc => grc_cs%beg_totsomc , & ! Input: [real(r8) (:)] (gC/m2) total column soil organic matter carbon
beg_totlitc => grc_cs%beg_totlitc , & ! Input: [real(r8) (:)] (gC/m2) total column litter carbon
beg_totprodc => grc_cs%beg_totprodc , & ! Input: [real(r8) (:)] (gC/m2) total column wood product carbon
beg_ctrunc => grc_cs%beg_ctrunc , & ! Input: [real(r8) (:)] (gC/m2) total column truncation carbon sink
beg_cropseedc_deficit => grc_cs%beg_cropseedc_deficit & ! Input: [real(r8) (:)] (gC/m2) column carbon pool for seeding new growth
)

! c2g states
Expand Down Expand Up @@ -946,13 +949,21 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
call c2g(bounds, col_som_c_leached(bounds%begc:bounds%endc), grc_som_c_leached(bounds%begg:bounds%endg), &
c2l_scale_type = 'unity', l2g_scale_type = 'unity')
call c2g(bounds, col_som_c_yield(bounds%begc:bounds%endc), grc_som_c_yield(bounds%begg:bounds%endg), &
c2l_scale_type = 'unity', l2g_scale_type = 'unity')
c2l_scale_type = 'unity', l2g_scale_type = 'unity')

if (use_fates) then
call c2g(bounds, col_cf%litfall(bounds%begc:bounds%endc), grc_cinputs(bounds%begg:bounds%endg), &
c2l_scale_type = 'unity', l2g_scale_type = 'unity')
end if

dt = real( get_step_size(), r8 )
nstep = get_nstep()

do g = bounds%begg, bounds%endg
grc_cinputs(g) = grc_gpp(g) + grc_dwt_seedc_to_leaf(g) + grc_dwt_seedc_to_deadstem(g)

if (.not. use_fates) then
grc_cinputs(g) = grc_gpp(g) + grc_dwt_seedc_to_leaf(g) + grc_dwt_seedc_to_deadstem(g)
end if

grc_coutputs(g) = grc_er(g) + grc_fire_closs(g) + grc_hrv_xsmrpool_to_atm(g) + &
grc_prod1c_loss(g) + grc_prod10c_loss(g) + grc_prod100c_loss(g) - grc_som_c_leached(g) + &
Expand All @@ -964,7 +975,7 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)

grc_errcb(g) = (grc_cinputs(g) - grc_coutputs(g))*dt - (end_totc(g) - beg_totc(g))

if (grc_errcb(g) > error_tol .and. nstep > 1) then
if (abs(grc_errcb(g)) > error_tol .and. nstep > 1) then
write(iulog,*)'grid cbalance error = ', grc_errcb(g), g
write(iulog,*)'Latdeg,Londeg = ', grc_pp%latdeg(g), grc_pp%londeg(g)
write(iulog,*)'input = ', grc_cinputs(g)*dt
Expand All @@ -973,6 +984,11 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
write(iulog,*)'fire = ', grc_fire_closs(g)*dt
write(iulog,*)'hrv_to_atm = ', grc_hrv_xsmrpool_to_atm(g)*dt
write(iulog,*)'leach = ', grc_som_c_leached(g)*dt
write(iulog,*)'prod1c_loss = ', grc_prod1c_loss(g)*dt
write(iulog,*)'prod10c_loss = ', grc_prod10c_loss(g)*dt
write(iulog,*)'prod100c_loss = ', grc_prod100c_loss(g)*dt
write(iulog,*)'som_c_leached = ', grc_som_c_leached(g)*dt
write(iulog,*)'dwt_conv_cflux = ', grc_dwt_conv_cflux(g)*dt

if (ero_ccycle) then
write(iulog,*)'erosion = ',grc_som_c_yield(g)*dt
Expand All @@ -982,6 +998,14 @@ subroutine GridCBalanceCheck(bounds, col_cs, col_cf, grc_cs, grc_cf)
write(iulog,*)'endcb = ', end_totc(g)
write(iulog,*)'delta store = ', end_totc(g) - beg_totc(g)

write(iulog,*)'Delta totpftc = ', end_totpftc(g) - beg_totpftc(g)
write(iulog,*)'Delta cwdc = ', end_cwdc(g) - beg_cwdc(g)
write(iulog,*)'Delta totlitc = ', end_totlitc(g) - beg_totlitc(g)
write(iulog,*)'Delta totsomc = ', end_totsomc(g) - beg_totsomc(g)
write(iulog,*)'Delta totprodc = ', end_totprodc(g) - beg_totprodc(g)
write(iulog,*)'Delta ctrunc = ', end_ctrunc(g) - beg_ctrunc(g)
write(iulog,*)'Delta crop deficit = ', end_cropseedc_deficit(g) - beg_cropseedc_deficit(g)

call endrun(msg=errMsg(__FILE__, __LINE__))
end if
end do
Expand Down
11 changes: 10 additions & 1 deletion components/elm/src/biogeochem/PhenologyMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2559,8 +2559,17 @@ subroutine CropPlantDate (num_soilp, filter_soilp, num_pcropp, filter_pcropp, &
call calculate_eto(t_ref2m(p), netrad, eflx_soil_grnd(p), forc_pbot(t), forc_rh(t), forc_wind(t), dt, ETout)
! monthly ETo
ETo(p,kmo) = ETo(p,kmo) + ETout

! calculate the P:PET for each month
p2ETo(p,kmo) = xp(p,kmo)/ETo(p,kmo)
if ( abs(ETo(p,kmo)) > 0._r8) then
p2ETo(p,kmo) = xp(p,kmo)/ETo(p,kmo)
else ! P:PET is undefined.
! Setting to a fill value ( 'spval' ) would
! require nested if statements due to
! the weighting of previous years (i.e., p2ETo and prev_p2ETo_bar )
! So, set to zero for simplicity.
p2ETo(p,kmo) = 0._r8
end if

if (nyrs_crop_active(p) == 0) then ! for the first year, use last years values
prev_xt_bar(p,kmo) = xt(p,kmo)
Expand Down
34 changes: 14 additions & 20 deletions components/elm/src/data_types/ColumnDataType.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7736,9 +7736,11 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column)
do fi = 1,num_column
i = filter_column(fi)
this%decomp_cpools_sourcesink(i,j,k) = value_column
this%decomp_cpools_transport_tendency(i,j,k) = value_column
end do
end do
end do

! pflotran
if(nstep_mod == 0 .or. is_first_restart_step()) then
do k = 1, ndecomp_pools
Expand Down Expand Up @@ -7776,15 +7778,17 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column)
this%er(i) = value_column
this%som_c_leached(i) = value_column
this%somc_yield(i) = value_column
this%somhr(i) = value_column ! REVISIT
this%lithr(i) = value_column ! REVISIT
this%somhr(i) = value_column
this%lithr(i) = value_column
this%hr(i) = value_column
this%cinputs(i) = value_column
this%coutputs(i) = value_column
this%cwdc_hr(i) = value_column
this%litterc_loss(i) = value_column

this%nee(i) = value_column
this%er(i) = value_column
this%som_c_leached(i) = value_column

! Zero p2c column fluxes
this%rr(i) = value_column
Expand All @@ -7804,36 +7808,26 @@ subroutine col_cf_setvalues ( this, num_column, filter_column, value_column)
this%somc_fire(i) = value_column
this%product_closs(i) = value_column
this%sr(i) = value_column
this%er(i) = value_column
this%litfire(i) = value_column
this%somfire(i) = value_column
this%totfire(i) = value_column
this%nep(i) = value_column
this%nbp(i) = value_column
this%fire_closs(i) = value_column
this%cwdc_loss(i) = value_column
this%som_c_leached(i) = value_column
this%somc_erode(i) = value_column
this%somc_deposit(i) = value_column
this%somc_yield(i) = value_column
enddo
end if
do k = 1, ndecomp_pools
do fi = 1,num_column
i = filter_column(fi)
this%decomp_cpools_leached(i,k) = value_column
this%decomp_cpools_erode(i,k) = value_column
this%decomp_cpools_deposit(i,k) = value_column
end do

do k = 1, ndecomp_pools
do fi = 1,num_column
i = filter_column(fi)
this%decomp_cpools_leached(i,k) = value_column
this%decomp_cpools_erode(i,k) = value_column
this%decomp_cpools_deposit(i,k) = value_column
end do
do k = 1, ndecomp_pools
do j = 1, nlevdecomp_full
do fi = 1,num_column
i = filter_column(fi)
this%decomp_cpools_transport_tendency(i,j,k) = value_column
end do
end do
end do
end do

end subroutine col_cf_setvalues

Expand Down

0 comments on commit 196a9b2

Please sign in to comment.