Skip to content

post5.4: Use the correct C13/C14 atm ratio in the code more #3632

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

Draft
ekluzek wants to merge 1 commit into
base: master
Choose a base branch
from
+27 −10
Draft

post5.4: Use the correct C13/C14 atm ratio in the code more #3632

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
16 changes: 11 additions & 5 deletions src/biogeochem/CNPhenologyMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ module CNPhenologyMod
use atm2lndType , only : atm2lnd_type
use CNVegMatrixMod , only : matrix_update_phc, matrix_update_phn
use CNVegMatrixMod , only : matrix_update_gmc, matrix_update_gmn
use CIsoAtmTimeseriesMod, only : C14BombSpike, nsectors_c14, C13TimeSeries
!
implicit none
private
Expand Down Expand Up @@ -128,6 +129,9 @@ module CNPhenologyMod
real(r8) :: hti ! cold hardening index threshold for vernalization
real(r8) :: tbase ! base temperature for vernalization

real(r8) :: rc14_atm(nsectors_c14)
real(r8) :: rc13_atm

integer, parameter :: NOT_Planted = 999 ! If not planted yet in year
integer, parameter :: NOT_Harvested = 999 ! If not harvested yet in year
integer, parameter :: inNH = 1 ! Northern Hemisphere
Expand Down Expand Up @@ -337,6 +341,7 @@ subroutine CNPhenology (bounds, num_soilc, filter_soilc, num_soilp, &
! !USES:
use clm_time_manager , only: is_first_step
use CNSharedParamsMod, only: use_fun
use clm_varctl , only : use_c13, use_c14
!
! !DESCRIPTION:
! Dynamic phenology routine for coupled carbon-nitrogen code (CN)
Expand Down Expand Up @@ -373,10 +378,12 @@ subroutine CNPhenology (bounds, num_soilc, filter_soilc, num_soilp, &
SHR_ASSERT_ALL_FL((ubound(leaf_prof_patch) == (/bounds%endp,nlevdecomp_full/)), sourcefile, __LINE__)
SHR_ASSERT_ALL_FL((ubound(froot_prof_patch) == (/bounds%endp,nlevdecomp_full/)), sourcefile, __LINE__)

if ( use_c13 ) call C13TimeSeries( rc13_atm )
if ( use_c14 ) call C14BombSpike( rc14_atm )

! each of the following phenology type routines includes a filter
! to operate only on the relevant patches


if ( phase == 1 ) then
call CNPhenologyClimate(num_soilp, filter_soilp, &
temperature_inst, cnveg_state_inst, crop_inst)
Expand Down Expand Up @@ -2014,7 +2021,6 @@ subroutine CropPhenology(num_pcropp, filter_pcropp , &
use clm_varcon , only : spval, secspday
use clm_varctl , only : use_fertilizer
use clm_varctl , only : use_c13, use_c14
use clm_varcon , only : c13ratio, c14ratio
use clm_varctl , only : use_cropcal_rx_swindows
!
! !ARGUMENTS:
Expand Down Expand Up @@ -2767,7 +2773,6 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, &
! !USES:
use clm_varctl , only : use_c13, use_c14
use clm_varctl , only : use_cropcal_rx_cultivar_gdds, adapt_cropcal_rx_cultivar_gdds
use clm_varcon , only : c13ratio, c14ratio
use clm_varpar , only : mxsowings
use pftconMod , only : ntmp_corn, nswheat, nwwheat, ntmp_soybean
use pftconMod , only : nirrig_tmp_corn, nirrig_swheat, nirrig_wwheat, nirrig_tmp_soybean
Expand Down Expand Up @@ -2864,15 +2869,16 @@ subroutine PlantCrop(p, leafcn_in, jday, kyr, do_plant_normal, &
c13_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * &
c13_cnveg_carbonstate_inst%totvegc_patch(p) / cnveg_carbonstate_inst%totvegc_patch(p)
else
c13_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * c13ratio
c13_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * rc13_atm
Copy link
Contributor

@wwieder wwieder Jan 28, 2026

Choose a reason for hiding this comment

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

This is only called if totvegC <0? I don't really understand how you can have leaf transfer C pool if totvegc is negative?

endif
endif
if (use_c14) then
if ( cnveg_carbonstate_inst%totvegc_patch(p) .gt. 0._r8) then
c14_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * &
c14_cnveg_carbonstate_inst%totvegc_patch(p) / cnveg_carbonstate_inst%totvegc_patch(p)
else
c14_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * c14ratio
! TODO: This should use the proper gridcell values rather than the global average
Copy link
Contributor

@wwieder wwieder Jan 28, 2026

Choose a reason for hiding this comment

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

As above, I'm not really clear when / why this would be called?

c14_cnveg_carbonstate_inst%leafc_xfer_patch(p) = leafc_xfer(p) * sum(rc14_atm(:))/real(nsectors_c14, r8)
endif
endif

Expand Down
8 changes: 6 additions & 2 deletions src/biogeochem/dynConsBiogeochemMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,10 +59,10 @@ subroutine dyn_cnbal_patch(bounds, &
use shr_const_mod , only : SHR_CONST_PDB
use landunit_varcon , only : istsoil, istcrop
use clm_varpar , only : nlevdecomp, i_litr_min, i_litr_max
use clm_varcon , only : c13ratio, c14ratio, c3_r2, c4_r2
use clm_time_manager , only : get_step_size_real
use dynPriorWeightsMod , only : prior_weights_type
use dynPatchStateUpdaterMod, only : patch_state_updater_type
use CIsoAtmTimeseriesMod, only : C13TimeSeries
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
Expand Down Expand Up @@ -126,6 +126,8 @@ subroutine dyn_cnbal_patch(bounds, &
real(r8), allocatable :: wood_product_c14flux(:) ! patch-level mass loss due to weight shift (expressed per unit GRIDCELL area)
real(r8), allocatable :: crop_product_c14flux(:) ! patch-level mass loss due to weight shift (expressed per unit GRIDCELL area)

real(r8) :: rc13_atm ! atmospheric C13/C12 ratio

logical :: patch_initiating(bounds%begp:bounds%endp)

! amounts to add to growing patches
Expand Down Expand Up @@ -308,6 +310,8 @@ subroutine dyn_cnbal_patch(bounds, &
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
endif

if ( use_c13 ) call C13TimeSeries( rc13_atm )

! Get time step
dt = get_step_size_real()
Expand Down Expand Up @@ -423,7 +427,7 @@ subroutine dyn_cnbal_patch(bounds, &
cnveg_nitrogenflux_inst%plant_nalloc_patch(p) = 0._r8

if ( use_c13 ) then
c13_cnveg_carbonflux_inst%xsmrpool_c13ratio_patch(p) = c13ratio
c13_cnveg_carbonflux_inst%xsmrpool_c13ratio_patch(p) = rc13_atm
end if

call photosyns_inst%NewPatchinit(p)
Expand Down
13 changes: 10 additions & 3 deletions src/soilbiogeochem/CNSoilMatrixMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ module CNSoilMatrixMod
use clm_varpar , only : ndecomp_pools, nlevdecomp, ndecomp_pools_vr !number of biogeochemically active soil layers
use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_cascade_outtransitions
use clm_varpar , only : i_cwd
use clm_varcon , only : dzsoi_decomp,zsoi,secspday,c3_r2,c14ratio
use clm_varcon , only : dzsoi_decomp,zsoi,secspday,c3_r2
use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, use_soil_matrixcn
use CNVegCarbonFluxType , only : cnveg_carbonflux_type
use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type
Expand All @@ -42,6 +42,7 @@ module CNSoilMatrixMod
use perf_mod , only : t_startf, t_stopf
use SparseMatrixMultiplyMod , only : sparse_matrix_type, diag_matrix_type, vector_type
use MatrixMod , only : inverse
use CIsoAtmTimeseriesMod, only : C14BombSpike, nsectors_c14
!
implicit none
private
Expand Down Expand Up @@ -140,6 +141,8 @@ subroutine CNSoilMatrix(bounds,num_soilc, filter_soilc, num_actfirec, filter_act
logical,save :: init_readyAsoilc = .False.
logical,save :: init_readyAsoiln = .False.
logical isbegofyear
real(r8) :: rc14_atm_sectors(nsectors_c14)
real(r8) :: rc14_atm

!-----------------------------------------------------------------------
begc = bounds%begc; endc = bounds%endc
Expand Down Expand Up @@ -224,6 +227,9 @@ subroutine CNSoilMatrix(bounds,num_soilc, filter_soilc, num_actfirec, filter_act
list_V_AKV => decomp_cascade_con%list_V_AKV &!In/Output:[Integer(:)] Saves mapping indices from V to (A*K+V) in the addition subroutine SPMP_AB
)

if ( use_c14 ) call C14BombSpike( rc14_atm_sectors )
rc14_atm = sum(rc14_atm_sectors(:)) / real(nsectors_c14,r8)

! set time steps
call t_startf('CN Soil matrix-init. matrix')
dt = real( get_step_size(), r8 )
Expand Down Expand Up @@ -362,13 +368,14 @@ subroutine CNSoilMatrix(bounds,num_soilc, filter_soilc, num_actfirec, filter_act
cs_soil%decomp0_cpools_vr_col = epsi
end where
if(use_c13)then
! Does the c3_r2 here need to change?
Copy link
Contributor

@wwieder wwieder Jan 28, 2026

Choose a reason for hiding this comment

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

This seems confusing, since c3_r2 is the "isotope ratio (13c/[12c+13c]) for C3 photosynthesis". This would ignore contribution from the 13C isotopic signature of C4 photosynthesis (c4_r2) to soil C pools?

Copy link
Contributor

@wwieder wwieder Jan 28, 2026

Choose a reason for hiding this comment

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

That said, this is also only called when decomp0_cpools_vr_col is a small number. Is this just for numerica stability?

where(cs13_soil%decomp0_cpools_vr_col .lt. epsi*c3_r2)
cs13_soil%decomp0_cpools_vr_col = epsi*c3_r2
end where
end if
if(use_c14)then
where(cs14_soil%decomp0_cpools_vr_col .lt. epsi*c14ratio)
cs14_soil%decomp0_cpools_vr_col = epsi*c14ratio
where(cs14_soil%decomp0_cpools_vr_col .lt. epsi*rc14_atm)
cs14_soil%decomp0_cpools_vr_col = epsi*rc14_atm
end where
end if
where(ns_soil%decomp0_npools_vr_col .lt. epsi)
Expand Down