Skip to content

Update Carbon isotope data to CMIP7 #3465

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
wants to merge 14 commits into
base: alpha-ctsm5.4.CMIP7
Choose a base branch
from
Draft

Update Carbon isotope data to CMIP7 #3465

Show file tree
Hide file tree
Changes from 13 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
fde0c08
Update nsectors to 4 so it will be hardcoded to work with just the ne…
ekluzek Aug 30, 2025
089c90e
Put in place the raw CMIP17 C isotope datasets from Imperial College,…
ekluzek Sep 1, 2025
bf5af0d
Merge branch 'ctsm5.4_update5.3.071' into CMIP7_carbon_isotopes
ekluzek Sep 2, 2025
c763c61
Hardcode the changes to get it to work with the new CMIP7 dataset
ekluzek Sep 2, 2025
9b1796d
Merge remote-tracking branch 'escomp/alpha-ctsm5.4.CMIP7' into CMIP7_…
ekluzek Sep 11, 2025
dab1a82
Reverse the loop for the CMIP7 data using beg and end indices and ste…
ekluzek Sep 11, 2025
363fa7f
Add RC14_CANAIR as an output history variable, and change the hardcod…
ekluzek Sep 17, 2025
e08b281
Set FC13/RC14 to spval so that zero's won't get averaged in
ekluzek Sep 18, 2025
8ed4536
Merge tag 'alpha-ctsm5.4.CMIP7.11.ctsm5.3.075' into CMIP7_carbon_isot…
ekluzek Sep 22, 2025
a0292c6
Add some logging about file times to see what's going on
ekluzek Sep 23, 2025
f64cf0d
Change the units of time so it's CE year to make it easier to underst…
ekluzek Sep 23, 2025
8117aba
Fix handling of sector 3
ekluzek Sep 23, 2025
b7786b5
Fix the latitude sectors for C14 so that it properly goes from South …
ekluzek Oct 10, 2025
30cef5e
Remove the incorrect comment noticed in the review
ekluzek Oct 10, 2025
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
4 changes: 2 additions & 2 deletions bld/namelist_files/namelist_defaults_ctsm.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,15 +122,15 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
<use_c14_bombspike phys="clm6_0" bgc_mode="bgc" >.true.</use_c14_bombspike>

<!-- Carbon isotope concentration files -->
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="hist" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_1850-2015_yearly_v2.0_c190528.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="hist" >lnd/clm2/isotopes/rawdata_copies/delta13co2_input4MIPs_atmosphericState_C4MIP_ImperialCollege-3-0_gm_1700-2023.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP1-1.9" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP119_1850-2100_yearly_c181209.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP1-2.6" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP126_1850-2100_yearly_c181209.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP2-4.5" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP245_1850-2100_yearly_c181209.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP3-7.0" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP3B_1850-2100_yearly_c181209.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP5-3.4" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP534os_1850-2100_yearly_c181209.nc</atm_c13_filename>
<atm_c13_filename use_c13=".true." use_c13_timeseries=".true." ssp_rcp="SSP5-8.5" >lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP5B_1850-2100_yearly_c181209.nc</atm_c13_filename>

<atm_c14_filename use_c14=".true." use_c14_bombspike =".true." ssp_rcp="hist" >lnd/clm2/isotopes/atm_delta_C14_CMIP6_3x1_global_1850-2015_yearly_v2.0_c190528.nc</atm_c14_filename>
<atm_c14_filename use_c14=".true." use_c14_bombspike =".true." ssp_rcp="hist" >lnd/clm2/isotopes/rawdata_copies/delta14co2_input4MIPs_atmosphericState_C4MIP_ImperialCollege-3-0_gz_1700-2023.nc</atm_c14_filename>
<atm_c14_filename use_c14=".true." use_c14_bombspike =".true." ssp_rcp="SSP1-1.9" >lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP119_3x1_global_1850-2100_yearly_c181209.nc</atm_c14_filename>
<atm_c14_filename use_c14=".true." use_c14_bombspike =".true." ssp_rcp="SSP1-2.6" >lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP126_3x1_global_1850-2100_yearly_c181209.nc</atm_c14_filename>
<atm_c14_filename use_c14=".true." use_c14_bombspike =".true." ssp_rcp="SSP2-4.5" >lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP245_3x1_global_1850-2100_yearly_c181209.nc</atm_c14_filename>
Expand Down
62 changes: 48 additions & 14 deletions src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ module CIsoAtmTimeseriesMod
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use clm_time_manager , only : get_curr_date, get_curr_yearfrac
use clm_time_manager , only : get_curr_date, get_curr_yearfrac, get_average_days_per_year
use clm_varcon , only : c14ratio, secspday
use shr_const_mod , only : SHR_CONST_PDB ! Ratio of C13/C12
use clm_varctl , only : iulog
Expand All @@ -27,7 +27,9 @@ module CIsoAtmTimeseriesMod
character(len=256) , public :: atm_c14_filename = ' ' ! file name of C14 input data
logical , public :: use_c13_timeseries = .false. ! do we use time-varying atmospheric C13?
character(len=256) , public :: atm_c13_filename = ' ' ! file name of C13 input data
integer, parameter , public :: nsectors_c14 = 3 ! Number of latitude sectors the C14 data has
integer, parameter , public :: nsectors_c14 = 4 ! Number of latitude sectors the C14 data has
integer, parameter , public :: beg_sector_c14 = nsectors_c14! The starting latitude sector for Northernmost latitude
integer, parameter , public :: end_sector_c14 = 1 ! The ending latitude sector for Southernnmost latitude

!
! !PRIVATE MEMBER FUNCTIONS:
Expand All @@ -38,7 +40,9 @@ module CIsoAtmTimeseriesMod
real(r8), allocatable, private :: atm_delta_c14(:,:) ! Delta C14 data
real(r8), allocatable, private :: atm_c13file_time(:) ! time for C13 data
real(r8), allocatable, private :: atm_delta_c13(:) ! Delta C13 data
real(r8), parameter :: time_axis_offset = 1850.0_r8 ! Offset in years of time on file
real(r8), parameter :: time_axis_offset = 1700.0_r8 ! Offset in years of time on file
real(r8), private :: c13_previous_time_idx = -1 ! Time index for C13 data on previous call
real(r8), private :: c14_previous_time_idx = -1 ! Time index for C14 data on previous call

character(len=*), parameter, private :: sourcefile = &
__FILE__
Expand All @@ -65,6 +69,7 @@ subroutine C14BombSpike( rc14_atm )
integer :: ntim_atm_ts ! Number of times on file
real(r8) :: twt_1, twt_2 ! weighting fractions for interpolating
integer :: l ! Loop index of sectors
integer :: sector_step = 1 ! Value to step the sector loop by (only 1 or -1)
!-----------------------------------------------------------------------

! get current date
Expand All @@ -75,13 +80,20 @@ subroutine C14BombSpike( rc14_atm )
ntim_atm_ts = size(atm_c14file_time)
ind_below = 0
do nt = 1, ntim_atm_ts
if ((dateyear - time_axis_offset) >= atm_c14file_time(nt) ) then
if ( dateyear >= atm_c14file_time(nt) ) then
ind_below = ind_below+1
endif
end do
! Save the previous time index and write out what index is being used if it changes
if ( (masterproc) .and. (c14_previous_time_idx /= ind_below) )then
write(iulog,*) ' c14 time index updated to just before current time is = ', ind_below, ' file year = ', &
atm_c14file_time(ind_below), ' current year = ', dateyear
end if
c14_previous_time_idx = ind_below

! loop over lat bands to pass all three to photosynthesis
do l = 1,nsectors_c14
if ( beg_sector_c14 > end_sector_c14 ) sector_step = -1
do l = beg_sector_c14, end_sector_c14, sector_step
! interpolate between nearest two points in atm c14 timeseries
if (ind_below .eq. 0 ) then
delc14o2_atm(l) = atm_delta_c14(l,1)
Expand Down Expand Up @@ -119,7 +131,7 @@ subroutine C14_init_BombSpike()
integer :: nsec ! number of input data sectors
integer :: t ! time index
logical :: readvar ! if variable read or not
character(len=*), parameter :: vname = 'Delta14co2_in_air' ! Variable name on file
character(len=*), parameter :: vname = 'Delta14co2' ! Variable name on file
!-----------------------------------------------------------------------

call getfil(atm_c14_filename, locfn, 0)
Expand All @@ -132,7 +144,7 @@ subroutine C14_init_BombSpike()
call ncd_pio_openfile (ncid, trim(locfn), 0)

call ncd_inqdlen(ncid,dimid,ntim,'time')
call ncd_inqdlen(ncid,dimid,nsec,'sector')
call ncd_inqdlen(ncid,dimid,nsec,'lat')
if ( nsec /= nsectors_c14 )then
call endrun(msg="ERROR: number of sectors on file not what's expected"//errMsg(sourcefile, __LINE__))
end if
Expand All @@ -157,13 +169,21 @@ subroutine C14_init_BombSpike()
call check_units( ncid, vname, "Modern" )
call ncd_pio_closefile(ncid)

! check to make sure that time dimension is well behaved
! Change units and check to make sure that time dimension is well behaved
atm_c14file_time(1) = time_axis_offset + (atm_c14file_time(1) / get_average_days_per_year())
do t = 2, ntim
! Convert time in days to years
atm_c14file_time(t) = time_axis_offset + (atm_c14file_time(t) / get_average_days_per_year())
if ( atm_c14file_time(t) - atm_c14file_time(t-1) <= 0._r8 ) then
write(iulog, *) 'C14_init_BombSpike: error. time axis must be monotonically increasing'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif
end do
if ( masterproc )then
write(iulog,*) ' c14 isotope file number of time samples is = ', ntim
write(iulog,*) ' c14 isotope, first year = ', atm_c14file_time(1)
write(iulog,*) ' c14 isotope, end year = ', atm_c14file_time(ntim)
end if

end subroutine C14_init_BombSpike

Expand Down Expand Up @@ -196,10 +216,16 @@ subroutine C13TimeSeries( rc13_atm )
ntim_atm_ts = size(atm_c13file_time)
ind_below = 0
do nt = 1, ntim_atm_ts
if ((dateyear - time_axis_offset) >= atm_c13file_time(nt) ) then
if ( dateyear >= atm_c13file_time(nt) ) then
ind_below = ind_below+1
endif
end do
! Save the previous time index and write out what index is being used if it changes
if ( (masterproc) .and. (c13_previous_time_idx /= ind_below) )then
write(iulog,*) ' c13 time index updated to just before current time is = ', ind_below, ' file year = ', atm_c13file_time(ind_below), &
' current year = ', dateyear
end if
c13_previous_time_idx = ind_below

! interpolate between nearest two points in atm c13 timeseries
! cdknotes. for now and for simplicity, just use the northern hemisphere values (sector 1)
Expand Down Expand Up @@ -238,7 +264,7 @@ subroutine C13_init_TimeSeries()
integer :: ntim ! number of input data time samples
integer :: t ! Time index
logical :: readvar ! if variable read or not
character(len=*), parameter :: vname = 'delta13co2_in_air' ! Variable name on file
character(len=*), parameter :: vname = 'delta13co2' ! Variable name on file
!-----------------------------------------------------------------------

call getfil(atm_c13_filename, locfn, 0)
Expand Down Expand Up @@ -272,13 +298,21 @@ subroutine C13_init_TimeSeries()
call check_units( ncid, vname, "VPDB" )
call ncd_pio_closefile(ncid)

! check to make sure that time dimension is well behaved
! Change units and check to make sure that time dimension is well behaved
atm_c13file_time(1) = time_axis_offset + (atm_c13file_time(1)/get_average_days_per_year())
do t = 2, ntim
atm_c13file_time(t) = time_axis_offset + (atm_c13file_time(t)/get_average_days_per_year())
if ( atm_c13file_time(t) - atm_c13file_time(t-1) <= 0._r8 ) then
write(iulog, *) 'C13_init_TimeSeries: error. time axis must be monotonically increasing'
call endrun(msg=errMsg(sourcefile, __LINE__))
endif
end do
if ( masterproc )then
write(iulog,*) ' c13 isotope file number of time samples is = ', ntim
write(iulog,*) ' c13 isotope, first year = ', atm_c13file_time(1)
write(iulog,*) ' c13 isotope, end year = ', atm_c13file_time(ntim)
end if


end subroutine C13_init_TimeSeries

Expand All @@ -304,14 +338,14 @@ subroutine check_units( ncid, vname, relativeto )

call ncd_inqvid( ncid, 'time', varid, vardesc )
call ncd_getatt( ncid, varid, "units", units )
write(t_units_expected,'("years since", I5, "-01-01 0:0:0.0")' ) nint(time_axis_offset)
write(t_units_expected,'("days since", I5, "-01-01 00:00:00")' ) nint(time_axis_offset)
if ( trim(units) /= t_units_expected )then
call endrun(msg="ERROR: time units on file are NOT what's expected"//errMsg(sourcefile, __LINE__))
end if
call ncd_inqvid( ncid, vname, varid, vardesc )
call ncd_getatt( ncid, varid, "units", units )
if ( trim(units) /= "per mil, relative to "//relativeto )then
call endrun(msg="ERROR: units on file for "//vname//" are NOT what's expected"//errMsg(sourcefile, __LINE__))
if ( trim(units) /= "1" )then
call endrun(msg="ERROR: units on file for "//vname//" are NOT what's expected", file=sourcefile, line=__LINE__ )
end if
end subroutine check_units

Expand Down
42 changes: 35 additions & 7 deletions src/biogeophys/PhotosynthesisMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,8 @@ module PhotosynthesisMod
real(r8), pointer, public :: rc13_psnsun_patch (:) ! patch C13O2/C12O2 in sunlit canopy psn flux
real(r8), pointer, public :: rc13_psnsha_patch (:) ! patch C13O2/C12O2 in shaded canopy psn flux

real(r8), pointer, public :: rc14_canair_patch (:) ! patch C14O2/C12O2 in canopy air

real(r8), pointer, public :: psnsun_patch (:) ! patch sunlit leaf photosynthesis (umol CO2/m**2/s)
real(r8), pointer, public :: psnsha_patch (:) ! patch shaded leaf photosynthesis (umol CO2/m**2/s)
real(r8), pointer, public :: c13_psnsun_patch (:) ! patch c13 sunlit leaf photosynthesis (umol 13CO2/m**2/s)
Expand Down Expand Up @@ -346,6 +348,7 @@ subroutine InitAllocate(this, bounds)
allocate(this%rc13_canair_patch (begp:endp)) ; this%rc13_canair_patch (:) = nan
allocate(this%rc13_psnsun_patch (begp:endp)) ; this%rc13_psnsun_patch (:) = nan
allocate(this%rc13_psnsha_patch (begp:endp)) ; this%rc13_psnsha_patch (:) = nan
allocate(this%rc14_canair_patch (begp:endp)) ; this%rc14_canair_patch (:) = nan

allocate(this%cisun_z_patch (begp:endp,1:nlevcan)) ; this%cisun_z_patch (:,:) = nan
allocate(this%cisha_z_patch (begp:endp,1:nlevcan)) ; this%cisha_z_patch (:,:) = nan
Expand Down Expand Up @@ -449,6 +452,7 @@ subroutine Clean(this)
deallocate(this%rc13_canair_patch )
deallocate(this%rc13_psnsun_patch )
deallocate(this%rc13_psnsha_patch )
deallocate(this%rc14_canair_patch )

deallocate(this%cisun_z_patch )
deallocate(this%cisha_z_patch )
Expand Down Expand Up @@ -579,7 +583,7 @@ subroutine InitHistory(this, bounds)
this%rc13_canair_patch(begp:endp) = spval
call hist_addfld1d (fname='RC13_CANAIR', units='proportion', &
avgflag='A', long_name='C13/C(12+13) for canopy air', &
ptr_patch=this%rc13_canair_patch, default='inactive')
ptr_patch=this%rc13_canair_patch, set_spec=spval, default='inactive')

this%rc13_psnsun_patch(begp:endp) = spval
call hist_addfld1d (fname='RC13_PSNSUN', units='proportion', &
Expand All @@ -592,6 +596,13 @@ subroutine InitHistory(this, bounds)
ptr_patch=this%rc13_psnsha_patch, default='inactive')
endif

if ( use_c14 ) then
this%rc14_canair_patch(begp:endp) = spval
call hist_addfld1d (fname='RC14_CANAIR', units='proportion', &
avgflag='A', long_name='C14/C(12+13) for canopy air', &
ptr_patch=this%rc14_canair_patch, set_spec=spval, default='inactive')
end if

! Canopy physiology

if ( use_c13 ) then
Expand Down Expand Up @@ -1034,6 +1045,11 @@ subroutine Restart(this, bounds, ncid, flag)
dim1name='pft', long_name='', units='', &
interpinic_flag='interp', readvar=readvar, data=this%rc13_psnsha_patch)
endif
if ( use_c14 ) then
call restartvar(ncid=ncid, flag=flag, varname='rc14_canair', xtype=ncd_double, &
dim1name='pft', long_name='', units='', &
interpinic_flag='interp', readvar=readvar, data=this%rc14_canair_patch)
end if

call restartvar(ncid=ncid, flag=flag, varname='GSSUN', xtype=ncd_double, &
dim1name='pft', dim2name='levcan', switchdim=.true., &
Expand Down Expand Up @@ -1173,10 +1189,13 @@ subroutine TimeStepInit (this, bounds)
.or. lun%itype(l) == istice &
.or. lun%itype(l) == istwet) then
if (use_c13) then
this%rc13_canair_patch(p) = 0._r8
this%rc13_canair_patch(p) = spval
this%rc13_psnsun_patch(p) = 0._r8
this%rc13_psnsha_patch(p) = 0._r8
end if
if (use_c14) then
this%rc14_canair_patch(p) = spval
end if
end if
end do

Expand All @@ -1197,11 +1216,15 @@ subroutine NewPatchInit (this, p)
if ( use_c13 ) then
this%alphapsnsun_patch(p) = 0._r8
this%alphapsnsha_patch(p) = 0._r8
this%rc13_canair_patch(p) = 0._r8
this%rc13_canair_patch(p) = spval
this%rc13_psnsun_patch(p) = 0._r8
this%rc13_psnsha_patch(p) = 0._r8
endif

if ( use_c14 ) then
this%rc14_canair_patch(p) = spval
end if

this%psnsun_patch(p) = 0._r8
this%psnsha_patch(p) = 0._r8

Expand Down Expand Up @@ -2072,6 +2095,7 @@ subroutine PhotosynthesisTotal (fn, filterp, &
rc13_canair => photosyns_inst%rc13_canair_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in canopy air
rc13_psnsun => photosyns_inst%rc13_psnsun_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in sunlit canopy psn flux
rc13_psnsha => photosyns_inst%rc13_psnsha_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in shaded canopy psn flux
rc14_canair => photosyns_inst%rc14_canair_patch , & ! Output: [real(r8) (:) ] C1342/C12O2 in canopy air
alphapsnsun => photosyns_inst%alphapsnsun_patch , & ! Output: [real(r8) (:) ] fractionation factor in sunlit canopy psn flux
alphapsnsha => photosyns_inst%alphapsnsha_patch , & ! Output: [real(r8) (:) ] fractionation factor in shaded canopy psn flux
psnsun_wc => photosyns_inst%psnsun_wc_patch , & ! Output: [real(r8) (:) ] Rubsico-limited sunlit leaf photosynthesis (umol CO2 /m**2/ s)
Expand Down Expand Up @@ -2134,14 +2158,18 @@ subroutine PhotosynthesisTotal (fn, filterp, &
if ( use_c14 ) then

! determine latitute sector for radiocarbon bomb spike inputs
if ( grc%latdeg(g) .ge. 30._r8 ) then
sector_c14 = 1
else if ( grc%latdeg(g) .ge. -30._r8 ) then
if ( grc%latdeg(g) >= 30._r8 ) then
sector_c14 = 4
else if ( grc%latdeg(g) >= 0._r8 ) then
sector_c14 = 3
else if ( grc%latdeg(g) >= -30._r8 ) then
sector_c14 = 2
else
sector_c14 = 3
sector_c14 = 1
endif

rc14_canair(p) = rc14_atm(sector_c14)

c14_psnsun(p) = rc14_atm(sector_c14) * psnsun(p)
c14_psnsha(p) = rc14_atm(sector_c14) * psnsha(p)
endif
Expand Down