diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml
index 80c96947b9..d1bec8e015 100644
--- a/bld/namelist_files/namelist_defaults_ctsm.xml
+++ b/bld/namelist_files/namelist_defaults_ctsm.xml
@@ -122,7 +122,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
.true.
-lnd/clm2/isotopes/atm_delta_C13_CMIP6_1850-2015_yearly_v2.0_c190528.nc
+lnd/clm2/isotopes/rawdata_copies/delta13co2_input4MIPs_atmosphericState_C4MIP_ImperialCollege-3-0_gm_1700-2023.nc
lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP119_1850-2100_yearly_c181209.nc
lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP126_1850-2100_yearly_c181209.nc
lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP245_1850-2100_yearly_c181209.nc
@@ -130,7 +130,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP534os_1850-2100_yearly_c181209.nc
lnd/clm2/isotopes/atm_delta_C13_CMIP6_SSP5B_1850-2100_yearly_c181209.nc
-lnd/clm2/isotopes/atm_delta_C14_CMIP6_3x1_global_1850-2015_yearly_v2.0_c190528.nc
+lnd/clm2/isotopes/rawdata_copies/delta14co2_input4MIPs_atmosphericState_C4MIP_ImperialCollege-3-0_gz_1700-2023.nc
lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP119_3x1_global_1850-2100_yearly_c181209.nc
lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP126_3x1_global_1850-2100_yearly_c181209.nc
lnd/clm2/isotopes/atm_delta_C14_CMIP6_SSP245_3x1_global_1850-2100_yearly_c181209.nc
diff --git a/src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90 b/src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90
index 529a547e88..d321433776 100644
--- a/src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90
+++ b/src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90
@@ -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
@@ -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:
@@ -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__
@@ -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
@@ -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)
@@ -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)
@@ -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
@@ -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
@@ -196,13 +216,18 @@ 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)
if (ind_below .eq. 0 ) then
delc13o2_atm = atm_delta_c13(1)
elseif (ind_below .eq. ntim_atm_ts ) then
@@ -238,7 +263,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)
@@ -272,13 +297,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
@@ -304,14 +337,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
diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90
index a823bce548..0156e2ab72 100644
--- a/src/biogeophys/PhotosynthesisMod.F90
+++ b/src/biogeophys/PhotosynthesisMod.F90
@@ -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)
@@ -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
@@ -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 )
@@ -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', &
@@ -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
@@ -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., &
@@ -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
@@ -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
@@ -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)
@@ -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