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.

Sign up for GitHub

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

Fix C balance and FPE errors #6235

Merged
merged 4 commits into from
Jun 5, 2024
Merged

Fix C balance and FPE errors #6235

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
4850808
Fixed divide-by-zero in crop model.
peterdschwartz Feb 1, 2024
a111285
Adjusted GridCBalanceCheck when FATES is active to set gridcell carbo…
peterdschwartz Feb 12, 2024
e6097d2
Adjusted setting of p2ETo in PhenologyMod.F90
peterdschwartz May 23, 2024
7b9ce63
Moved fates_land_developer test to fates test suite until it can be f…
peterdschwartz Jun 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion cime_config/tests.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,6 @@
"fates_elm_developer" : {
"inherit" : ("fates_elm_debug"),
"tests" : (
"ERP_Ld15.ne4pg2_ne4pg2.IELMFATES.elm-fates_cold_allvars",
"ERS_Ld30.f45_f45.IELMFATES.elm-fates_satphen",
"ERS_Ld30.f45_g37.IELMFATES.elm-fates_cold_sizeagemort",
"SMS_Ld20.f45_f45.IELMFATES.elm-fates_eca",
Expand All @@ -437,6 +436,7 @@
"fates" : {
"inherit" : ("fates_long_tests", "fates_elm_developer"),
"tests" : (
"ERP_Ld15.ne4pg2_ne4pg2.IELMFATES.elm-fates_cold_allvars",
"ERP_Ld3.f09_g16.IELMFATES.elm-fates_cold",
"ERP_D_Ld3.f19_g16.IELMFATES.elm-fates_cold",
"ERS_D_Ld3_PS.f09_g16.IELMFATES.elm-fates_cold",
Expand Down
106 changes: 65 additions & 41 deletions components/elm/src/biogeochem/EcosystemBalanceCheckMod.F90
View file
Open in desktop
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
Open in desktop
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
Comment on lines +2564 to +2572
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@peterdschwartz This looks good to me. Thanks for the comment explaining the logic.


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
Open in desktop
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