[WIP]: Start bringing in carbon isotopes in streams

src/biogeochem/AtmCarbonIsotopeStreamType.F90

@@ -0,0 +1,193 @@
+module AtmCarbonIsotopeStreamType
+  use shr_kind_mod , only : r8 => shr_kind_r8
+  use clm_varctl , only : iulog
+  use abortutils , only : endrun
+  use decompMod , only : bounds_type
+  use CTSMForce2DStreamBaseType, only : ctsm_force_2DStream_base_type
+  use dshr_methods_mod , only : dshr_fldbun_getfldptr
+  use ESMF, only : ESMF_LogFoundError, ESMF_LOGERR_PASSTHRU
+
+  implicit none
+  private
+
+  character(len=*), parameter :: varname_c13 = 'delta13co2_in_air'
+  type, public, extends(ctsm_force_2DStream_base_type) :: atm_delta_c13_stream_type
+     private
+     real(r8), allocatable :: atm_delta_c13(:) ! delta C13 data array
+  contains
+
+      ! Public Methods
+      procedure, public :: C13Init
+      procedure, public :: Init => C13Init
+      procedure, public :: C13Interp 
+      procedure, public :: Interp => C13Interp 
+      procedure, public :: C13ClassClean
+      procedure, public :: Clean => C13ClassClean
+      final :: C13TypeClean
+      ! Private methods
+      procedure, private :: C13InitAllocate
+
+  end type atm_delta_c13_stream_type
+
+  character(len=*), parameter :: varname_c14 = 'Delta14co2_in_air'
+  type, public, extends(ctsm_force_2DStream_base_type) :: atm_delta_c14_stream_type
+     private
+     real(r8), allocatable :: atm_delta_c14(:) ! delta c14 data array
+  contains
+
+      ! Public Methods
+      procedure, public :: C14Init
+      procedure, public :: Init => C14Init
+      procedure, public :: C14Interp 
+      procedure, public :: Interp => C14Interp 
+      procedure, public :: C14ClassClean
+      procedure, public :: Clean => C14ClassClean
+      final :: C14TypeClean
+      ! Private methods
+      procedure, private :: C14InitAllocate
+
+  end type atm_delta_c14_stream_type
+
+  character(len=*), parameter, private :: sourcefile = &
+  __FILE__
+
+  contains
+
+    subroutine C13Init( this, bounds, fldfilename, meshfile, mapalgo, tintalgo, taxmode, &
+                        year_first, year_last, model_year_align )
+         ! Uses:
+         ! Arguments:
+         class(atm_delta_c13_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+         character(*), intent(in) :: fldfilename ! stream data filename (full pathname) (single file)
+         character(*), intent(in) :: meshfile ! full pathname to stream mesh file (none for global data)
+         character(*), intent(in) :: mapalgo ! stream mesh -> model mesh mapping type
+         character(*), intent(in) :: tintalgo ! time interpolation algorithm
+         character(*), intent(in) :: taxMode ! time axis mode
+         integer, intent(in) :: year_first ! first year to use
+         integer, intent(in) :: year_last ! last  year to use
+         integer, intent(in) :: model_year_align ! align yearFirst with this model year
+
+         call this%InitBase( bounds, varnames = (/ varname_c13 /), fldfilename=fldfilename, meshfile=meshfile, &
+                             mapalgo=mapalgo, tintalgo=tintalgo, taxmode=taxmode, name=varname_c13, &
+                             year_first=year_first, year_last=year_last, model_year_align=model_year_align )
+         call this%C13InitAllocate( bounds )
+
+     end subroutine C13Init
+
+    subroutine C13InitAllocate( this, bounds )
+         use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+         class(atm_delta_c13_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+
+         integer :: begg, endg
+
+         begg = bounds%begg; endg = bounds%endg
+         allocate( this%atm_delta_c13( bounds%begg : bounds%endg ) ); this%atm_delta_c13 = nan
+    end subroutine C13InitAllocate
+
+    subroutine C13ClassClean( this )
+         class(atm_delta_c13_stream_type), intent(inout) :: this 
+
+         call C13TypeClean( this )
+
+    end subroutine C13ClassClean
+
+    subroutine C13TypeClean( this )
+         type(atm_delta_c13_stream_type), intent(inout) :: this 
+
+         deallocate( this%atm_delta_c13 )
+         call this%CleanBase()
+
+    end subroutine C13TypeClean
+
+    subroutine C13Interp( this, bounds )
+         class(atm_delta_c13_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+
+         ! Local Variables
+         integer :: g
+         real(r8), pointer :: dataptr1d(:)
+         integer :: rc ! error return code
+
+         ! Get pointer for stream data that is time and spatially interpolated to model time and grid
+         call dshr_fldbun_getFldPtr(this%sdat%pstrm(1)%fldbun_model, fldname=varname_c13, fldptr1=dataptr1d, rc=rc)
+         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=sourcefile)) then
+            call endrun( 'Error getting field pointer for '//varname_c13//' from stream data', file=sourcefile, line=__LINE__ )
+         end if
+
+         do g = bounds%begg, bounds%endg
+            this%atm_delta_c13(g) = dataptr1d(g)
+         end do
+     end subroutine C13Interp
+
+     subroutine C14Init( this, bounds, fldfilename, meshfile, mapalgo, tintalgo, taxmode, &
+                         year_first, year_last, model_year_align )
+         ! Uses:
+         ! Arguments:
+         class(atm_delta_c14_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+         character(*), intent(in) :: fldfilename ! stream data filename (full pathname) (single file)
+         character(*), intent(in) :: meshfile ! full pathname to stream mesh file (none for global data)
+         character(*), intent(in) :: mapalgo ! stream mesh -> model mesh mapping type
+         character(*), intent(in) :: tintalgo ! time interpolation algorithm
+         character(*), intent(in) :: taxMode ! time axis mode
+         integer, intent(in) :: year_first ! first year to use
+         integer, intent(in) :: year_last ! last  year to use
+         integer, intent(in) :: model_year_align ! align yearFirst with this model year
+
+         call this%InitBase( bounds, varnames = (/ varname_c14 /), fldfilename=fldfilename, meshfile=meshfile, &
+                             mapalgo=mapalgo, tintalgo=tintalgo, taxmode=taxmode, name=varname_c14, &
+                             year_first=year_first, year_last=year_last, model_year_align=model_year_align )
+         call this%C14InitAllocate( bounds )
+
+     end subroutine C14Init
+
+    subroutine C14InitAllocate( this, bounds )
+         use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+         class(atm_delta_c14_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+
+         integer :: begg, endg
+
+         begg = bounds%begg; endg = bounds%endg
+         allocate( this%atm_delta_c14( bounds%begg : bounds%endg ) ); this%atm_delta_c14 = nan
+    end subroutine C14InitAllocate
+
+    subroutine C14Interp( this, bounds )
+         class(atm_delta_c14_stream_type), intent(inout) :: this 
+         type(bounds_type), intent(in) :: bounds
+
+         ! Local Variables
+         integer :: g
+         real(r8), pointer :: dataptr1d(:)
+         integer :: rc ! error return code
+
+         ! Get pointer for stream data that is time and spatially interpolated to model time and grid
+         call dshr_fldbun_getFldPtr(this%sdat%pstrm(1)%fldbun_model, fldname=varname_c14, fldptr1=dataptr1d, rc=rc)
+         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=sourcefile)) then
+            call endrun( 'Error getting field pointer for '//varname_c14//' from stream data', file=sourcefile, line=__LINE__ )
+         end if
+
+         do g = bounds%begg, bounds%endg
+            this%atm_delta_c14(g) = dataptr1d(g)
+         end do
+    end subroutine C14Interp
+
+    subroutine C14ClassClean( this )
+         class(atm_delta_c14_stream_type), intent(inout) :: this 
+
+         call C14TypeClean( this )
+
+    end subroutine C14ClassClean
+
+
+    subroutine C14TypeClean( this )
+         type(atm_delta_c14_stream_type), intent(inout) :: this 
+
+         deallocate( this%atm_delta_c14 )
+         call this%CleanBase()
+
+    end subroutine C14TypeClean
+
+end module AtmCarbonIsotopeStreamType

src/biogeochem/CNCIsoAtmTimeSeriesReadMod.F90

@@ -4,14 +4,16 @@ module CIsoAtmTimeseriesMod
   ! Module for transient atmospheric boundary to the c13 and c14 codes
   !
   ! !USES:
-  use shr_kind_mod     , only : r8 => shr_kind_r8
+  use shr_kind_mod     , only : r8 => shr_kind_r8, CL => shr_kind_CL
   use clm_time_manager , only : get_curr_date, get_curr_yearfrac
   use clm_varcon       , only : c14ratio, secspday
   use shr_const_mod    , only : SHR_CONST_PDB                    ! Ratio of C13/C12
   use clm_varctl       , only : iulog
   use abortutils       , only : endrun
   use spmdMod          , only : masterproc
   use shr_log_mod      , only : errMsg => shr_log_errMsg
+  use AtmCarbonIsotopeStreamType, only : atm_delta_c13_stream_type, atm_delta_c14_stream_type
+  use decompMod        , only : bounds_type
   !
   implicit none
   private
@@ -33,13 +35,33 @@ module CIsoAtmTimeseriesMod
   ! !PRIVATE MEMBER FUNCTIONS:
   private:: check_units   ! Check the units of the data on the input file
 
+  type(atm_delta_c13_stream_type), private :: atm_c13_stream  ! Atmospheric C13 stream object
+  type(atm_delta_c14_stream_type), private :: atm_c14_stream  ! Atmospheric C14 stream object
+
   ! !PRIVATE TYPES:
   real(r8), allocatable, private :: atm_c14file_time(:)       ! time for C14 data
   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
 
+  ! Private data for the control namelist:
+  character(len=CL), private :: stream_fldfilename_atm_c14 = '/glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/isotopes/ctsmforc.Graven.atm_delta_C14_CMIP7_4x1_global_1700-2023_yearly_v3.0_c251013.nc'
+  character(len=CL), private :: stream_meshfile_atm_c14 = '/glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/isotopes/mesh_4x1_global_c20251013.nc'
+  character(len=CL), private :: stream_fldfilename_atm_c13 = '/glade/campaign/cesm/cesmdata/cseg/inputdata/lnd/clm2/isotopes/ctsmforc.Graven.atm_delta_C13_CMIP7_global_1700-2023_yearly_v3.0_c251013.nc'
+  integer, private :: stream_year_first_atm_c14 = 1850
+  integer, private :: stream_year_last_atm_c14 = 2023
+  integer, private :: stream_model_year_align_atm_c14 = 1850
+  integer, private :: stream_year_first_atm_c13 = 1850
+  integer, private :: stream_year_last_atm_c13 = 2023
+  integer, private :: stream_model_year_align_atm_c13 = 1850
+  character(len=CL), private :: stream_mapalgo_atm_c14 = 'nn'
+  character(len=CL), private :: stream_tintalgo_atm_c14 = 'linear'
+  character(len=CL), private :: stream_taxmode_atm_c14 = 'extend'
+  character(len=CL), private :: stream_mapalgo_atm_c13 = 'nn'
+  character(len=CL), private :: stream_tintalgo_atm_c13 = 'linear'
+  character(len=CL), private :: stream_taxmode_atm_c13 = 'extend'
+
   character(len=*), parameter, private :: sourcefile = &
        __FILE__
   !-----------------------------------------------------------------------
@@ -109,6 +131,7 @@ subroutine C14_init_BombSpike()
     ! !USES:
     use ncdio_pio   , only : ncd_pio_openfile, ncd_pio_closefile, file_desc_t, ncd_inqdlen, ncd_io
     use fileutils   , only : getfil
+    use decompMod   , only : get_proc_bounds
     !
     ! !LOCAL VARIABLES:
     implicit none
@@ -119,6 +142,7 @@ subroutine C14_init_BombSpike()
     integer :: nsec                       ! number of input data sectors
     integer :: t                          ! time index
     logical :: readvar                    ! if variable read or not
+    type(bounds_type) :: bounds_proc
     character(len=*), parameter :: vname = 'Delta14co2_in_air'  ! Variable name on file
     !-----------------------------------------------------------------------
 
@@ -165,6 +189,21 @@ subroutine C14_init_BombSpike()
        endif
     end do
 
+    ! Streams method
+    call get_proc_bounds( bounds_proc )
+    call atm_c14_stream%Init( bounds_proc, &
+        fldfilename=stream_fldfilename_atm_c14, &
+        !meshfile= stream_meshfile_atm_c14, &
+        meshfile= 'none', &
+        mapalgo=stream_mapalgo_atm_c14, &
+        tintalgo=stream_tintalgo_atm_c14, &
+        taxmode=stream_taxmode_atm_c14, &
+        year_first=stream_year_first_atm_c14, &
+        year_last=stream_year_last_atm_c14, &
+        model_year_align=stream_model_year_align_atm_c14 )
+    call atm_c14_stream%Advance( )
+    call atm_c14_stream%Interp( bounds_proc )
+
   end subroutine C14_init_BombSpike
 
 
@@ -229,6 +268,7 @@ subroutine C13_init_TimeSeries()
     ! !USES:
     use ncdio_pio   , only : ncd_pio_openfile, ncd_pio_closefile, file_desc_t, ncd_inqdlen, ncd_io
     use fileutils   , only : getfil
+    use decompMod   , only : get_proc_bounds
     !
     ! !LOCAL VARIABLES:
     implicit none
@@ -238,6 +278,7 @@ subroutine C13_init_TimeSeries()
     integer :: ntim                       ! number of input data time samples
     integer :: t                          ! Time index
     logical :: readvar                    ! if variable read or not
+    type(bounds_type) :: bounds_proc
     character(len=*), parameter :: vname = 'delta13co2_in_air'  ! Variable name on file
     !-----------------------------------------------------------------------
 
@@ -280,6 +321,20 @@ subroutine C13_init_TimeSeries()
        endif
     end do
 
+    ! Streams method
+    call get_proc_bounds( bounds_proc )
+    call atm_c13_stream%Init( bounds_proc, &
+        fldfilename=stream_fldfilename_atm_c13, &
+        meshfile= 'none', &
+        mapalgo=stream_mapalgo_atm_c13, &
+        tintalgo=stream_tintalgo_atm_c13, &
+        taxmode=stream_taxmode_atm_c13, &
+        year_first=stream_year_first_atm_c13, &
+        year_last=stream_year_last_atm_c13, &
+        model_year_align=stream_model_year_align_atm_c13 )
+    call atm_c13_stream%Advance( )
+    call atm_c13_stream%Interp( bounds_proc )
+
   end subroutine C13_init_TimeSeries
 
   !-----------------------------------------------------------------------