diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.F90
new file mode 100644
index 000000000..fbaf5a1d9
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.F90
@@ -0,0 +1,442 @@
+! ########################################################################################
+! 
+! Description: Interstitial CCPP suite to couple UFS physics to ccpp_suite_simulator.
+!
+! Contains:
+! - load_ccpp_suite_sim(): read and load data into type used by ccpp_suite_simulator.
+!      called once during model initialization
+! - GFS_ccpp_suite_sim_pre_run(): prepare GFS diagnostic physics tendencies for 
+!      ccpp_suite_simulator. 
+!
+! ########################################################################################
+module GFS_ccpp_suite_sim_pre
+  use machine, only: kind_phys
+  use module_ccpp_suite_simulator, only: base_physics_process
+  use netcdf
+  implicit none
+  public GFS_ccpp_suite_sim_pre_run, load_ccpp_suite_sim
+contains
+
+  ! ######################################################################################
+  !
+  ! SUBROUTINE GFS_ccpp_suite_sim_pre_run
+  !
+  ! ######################################################################################
+!! \section arg_table_GFS_ccpp_suite_sim_pre_run
+!! \htmlinclude GFS_ccpp_suite_sim_pre_run.html
+!! 
+  subroutine GFS_ccpp_suite_sim_pre_run(do_ccpp_suite_sim, dtend, ntqv, dtidx, dtp,      &
+       index_of_process_dcnv, index_of_process_longwave, index_of_process_shortwave,     &
+       index_of_process_scnv, index_of_process_orographic_gwd, index_of_process_pbl,     &
+       index_of_process_mp, index_of_temperature, index_of_x_wind, index_of_y_wind,      &
+       physics_process, iactive_T, iactive_u, iactive_v, iactive_q, active_phys_tend,    &
+       errmsg, errflg)
+
+    ! Inputs
+    logical, intent(in) :: do_ccpp_suite_sim
+    integer, intent(in) :: ntqv, index_of_process_dcnv, index_of_process_longwave,       &
+         index_of_process_shortwave, index_of_process_scnv,                              &
+         index_of_process_orographic_gwd, index_of_process_pbl, index_of_process_mp,     &
+         index_of_temperature, index_of_x_wind, index_of_y_wind
+    integer, intent(in), dimension(:,:) :: dtidx
+    real(kind_phys), intent(in) :: dtp
+    real(kind_phys), intent(in), dimension(:,:,:) :: dtend
+    type(base_physics_process),intent(in) :: physics_process(:)
+    integer,         intent(in) :: iactive_T, iactive_u, iactive_v, iactive_q
+
+    ! Outputs
+    real(kind_phys), intent(out) :: active_phys_tend(:,:,:)
+    character(len=*),intent(out) :: errmsg
+    integer,         intent(out) :: errflg
+
+    ! Locals
+    integer :: idtend, iactive
+
+    ! Initialize CCPP error handling variables
+    errmsg = ''
+    errflg = 0
+
+    if (.not. do_ccpp_suite_sim) return
+
+    ! Get tendency for "active" process.
+
+    ! ######################################################################################
+    ! DJS2023: For the UFS and SCM, the physics tendencies are stored in a multi-dimensional
+    ! array, CCPP standard_name = cumulative_change_of_state_variables.
+    ! These are not the instantaneous physics tendencies that are applied to the state by 
+    ! the physics suites. Not all suites output physics tendencies...
+    ! Rather these are intended for diagnostic puposes and are accumulated over some 
+    ! interval.
+    ! In the UFS/SCM this is controlled by the diagnostic bucket interval, namelist option 
+    ! "fhzero". For this to work, you need to clear the diagnostic buckets after each
+    ! physics timestep when running in the UFS/SCM.
+    !
+    ! In the SCM this is done by adding the following runtime options:
+    ! --n_itt_out 1 --n_itt_diag 1
+    !
+    ! ######################################################################################
+    if (physics_process(1)%active_name == "LWRAD") iactive = index_of_process_longwave
+    if (physics_process(1)%active_name == "SWRAD") iactive = index_of_process_shortwave
+    if (physics_process(1)%active_name == "PBL")   iactive = index_of_process_pbl
+    if (physics_process(1)%active_name == "GWD")   iactive = index_of_process_orographic_gwd
+    if (physics_process(1)%active_name == "SCNV")  iactive = index_of_process_scnv
+    if (physics_process(1)%active_name == "DCNV")  iactive = index_of_process_dcnv
+    if (physics_process(1)%active_name == "cldMP") iactive = index_of_process_mp
+
+    ! Heat
+    idtend = dtidx(index_of_temperature,iactive)
+    if (idtend >= 1) then
+       active_phys_tend(:,:,iactive_T) = dtend(:,:,idtend)/dtp
+    endif
+
+    ! u-wind
+    idtend = dtidx(index_of_x_wind,iactive)
+    if (idtend >= 1) then
+       active_phys_tend(:,:,iactive_u) = dtend(:,:,idtend)/dtp
+    endif
+
+    ! v-wind
+    idtend = dtidx(index_of_y_wind,iactive)
+    if (idtend >= 1) then
+       active_phys_tend(:,:,iactive_v) = dtend(:,:,idtend)/dtp
+    endif
+
+    ! Moisture
+    idtend = dtidx(100+ntqv,iactive)
+    if (idtend >= 1) then
+       active_phys_tend(:,:,iactive_q) = dtend(:,:,idtend)/dtp
+    endif
+
+  end subroutine GFS_ccpp_suite_sim_pre_run
+
+  ! ######################################################################################
+  subroutine load_ccpp_suite_sim(nlunit, nml_file, physics_process, iactive_T,           &
+       iactive_u, iactive_v, iactive_q, errmsg, errflg)
+
+    ! Inputs
+    integer,          intent (in) :: nlunit
+    character(len=*), intent (in) :: nml_file
+
+    ! Outputs
+    type(base_physics_process),intent(inout),allocatable :: physics_process(:)
+    integer, intent(inout)        :: iactive_T, iactive_u, iactive_v, iactive_q
+    integer, intent(out)          :: errflg
+    character(len=256), intent(out) :: errmsg
+
+    ! Local variables
+    integer :: ncid, dimID, varID, status, ios, iprc, nlev_data, ntime_data
+    character(len=256) :: suite_sim_file
+    logical :: exists, do_ccpp_suite_sim
+    integer :: nprc_sim
+
+    ! For each process there is a corresponding namelist entry, which is constructed as 
+    ! follows:
+    ! {use_suite_sim[0(no)/1(yes)], time_split[0(no)/1(yes)], order[1:nPhysProcess]}
+    integer, dimension(3) ::    &
+         prc_LWRAD_cfg = (/0,0,0/), &
+         prc_SWRAD_cfg = (/0,0,0/), &
+         prc_PBL_cfg   = (/0,0,0/), &
+         prc_GWD_cfg   = (/0,0,0/), &
+         prc_SCNV_cfg  = (/0,0,0/), &
+         prc_DCNV_cfg  = (/0,0,0/), &
+         prc_cldMP_cfg = (/0,0,0/)
+
+    ! Namelist
+    namelist / ccpp_suite_sim_nml / do_ccpp_suite_sim, suite_sim_file, nprc_sim,         &
+         prc_LWRAD_cfg, prc_SWRAD_cfg, prc_PBL_cfg, prc_GWD_cfg, prc_SCNV_cfg,           &
+         prc_DCNV_cfg, prc_cldMP_cfg
+
+    errmsg = ''
+    errflg = 0
+
+    ! Read in namelist
+    inquire (file = trim (nml_file), exist = exists)
+    if (.not. exists) then
+        errmsg = 'CCPP suite simulator namelist file: '//trim(nml_file)//' does not exist.'
+        errflg = 1
+        return
+    else
+        open (unit = nlunit, file = nml_file, action = 'read', status = 'old', iostat = ios)
+    endif
+    rewind (nlunit)
+    read (nlunit, nml = ccpp_suite_sim_nml, iostat=status)
+    close (nlunit)
+
+    ! Only proceed if suite simulator requested.
+    if (prc_SWRAD_cfg(1)  == 1 .or. prc_LWRAD_cfg(1) == 1 .or. prc_PBL_cfg(1)  == 1 .or. &
+         prc_GWD_cfg(1)   == 1 .or. prc_SCNV_cfg(1)  == 1 .or. prc_DCNV_cfg(1) == 1 .or. &
+         prc_cldMP_cfg(1) == 1 ) then
+    else
+       return
+    endif
+
+    ! Check that input data file exists.
+    inquire (file = trim (suite_sim_file), exist = exists)
+    if (.not. exists) then
+       errmsg = 'CCPP suite simulator file: '//trim(suite_sim_file)//' does not exist'
+       errflg = 1
+       return
+    endif
+
+    !
+    ! Read data file...
+    !
+
+    ! Open file
+    status = nf90_open(trim(suite_sim_file), NF90_NOWRITE, ncid)
+    if (status /= nf90_noerr) then
+       errmsg = 'Error reading in CCPP suite simulator file: '//trim(suite_sim_file)
+       errflg = 1
+       return
+    endif
+
+    ! Metadata (dimensions)
+    status = nf90_inq_dimid(ncid, 'time', dimid)
+    if (status == nf90_noerr) then
+       status = nf90_inquire_dimension(ncid, dimid, len = ntime_data)
+    else
+       errmsg = 'CCPP suite simulator file: '//trim(suite_sim_file)//' does not contain [time] dimension'
+       errflg = 1
+       return
+    endif
+ 
+    status = nf90_inq_dimid(ncid, 'lev', dimid)
+    if (status == nf90_noerr) then
+       status = nf90_inquire_dimension(ncid, dimid, len = nlev_data)
+    else
+       errmsg = 'CCPP suite simulator file: '//trim(suite_sim_file)//' does not contain [lev] dimension'
+       errflg = 1
+       return
+    endif
+
+    ! Allocate space and read in data
+    allocate(physics_process(nprc_sim))
+    physics_process(1)%active_name    = ''
+    physics_process(1)%iactive_scheme = 0
+    physics_process(1)%active_tsp     = .false.
+    do iprc = 1,nprc_sim
+       allocate(physics_process(iprc)%tend1d%T(   nlev_data            ))
+       allocate(physics_process(iprc)%tend1d%u(   nlev_data            ))
+       allocate(physics_process(iprc)%tend1d%v(   nlev_data            ))
+       allocate(physics_process(iprc)%tend1d%q(   nlev_data            ))
+       allocate(physics_process(iprc)%tend2d%time(           ntime_data))
+       allocate(physics_process(iprc)%tend2d%T(   nlev_data, ntime_data))
+       allocate(physics_process(iprc)%tend2d%u(   nlev_data, ntime_data))
+       allocate(physics_process(iprc)%tend2d%v(   nlev_data, ntime_data))
+       allocate(physics_process(iprc)%tend2d%q(   nlev_data, ntime_data))
+
+       ! Temporal info
+       status = nf90_inq_varid(ncid, 'times', varID)
+       if (status == nf90_noerr) then
+          status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%time)
+       else
+          errmsg = 'SCM data tendency file: '//trim(suite_sim_file)//' does not contain times variable'
+          errflg = 1
+          return
+       endif
+
+       if (iprc == prc_SWRAD_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "SWRAD"
+          if (prc_SWRAD_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 1
+             iactive_T   = 1
+          endif
+          if (prc_SWRAD_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_swrad', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+       endif
+
+       if (iprc == prc_LWRAD_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "LWRAD"
+          if (prc_LWRAD_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 1
+             iactive_T   = 1
+          endif
+          if (prc_LWRAD_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_lwrad', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+       endif
+
+       if (iprc == prc_GWD_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "GWD"
+          if (prc_GWD_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 3
+             iactive_T   = 1
+             iactive_u   = 2
+             iactive_v   = 3
+          endif
+          if (prc_GWD_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_cgwd', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+          status = nf90_inq_varid(ncid, 'du_dt_cgwd', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%u)
+          status = nf90_inq_varid(ncid, 'dv_dt_cgwd', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%v)
+       endif
+
+       if (iprc == prc_PBL_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "PBL"
+          if (prc_PBL_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 4
+             iactive_T   = 1
+             iactive_u   = 2
+             iactive_v   = 3
+             iactive_q   = 4
+          endif
+          if (prc_PBL_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_pbl', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+          status = nf90_inq_varid(ncid, 'dq_dt_pbl', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%q)
+          status = nf90_inq_varid(ncid, 'du_dt_pbl', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%u)
+          status = nf90_inq_varid(ncid, 'dv_dt_pbl', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%v)
+       endif
+
+       if (iprc == prc_SCNV_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "SCNV"
+          if (prc_SCNV_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 4
+             iactive_T   = 1
+             iactive_u   = 2
+             iactive_v   = 3
+             iactive_q   = 4
+          endif
+          if (prc_SCNV_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_shalconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+          status = nf90_inq_varid(ncid, 'du_dt_shalconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%u)
+          status = nf90_inq_varid(ncid, 'dv_dt_shalconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%v)
+          status = nf90_inq_varid(ncid, 'dq_dt_shalconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%q)
+       endif
+
+       if (iprc == prc_DCNV_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "DCNV"
+          if (prc_DCNV_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 4
+             iactive_T   = 1
+             iactive_u   = 2
+             iactive_v   = 3
+             iactive_q   = 4
+          endif
+          if (prc_DCNV_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_deepconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+          status = nf90_inq_varid(ncid, 'du_dt_deepconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%u)
+          status = nf90_inq_varid(ncid, 'dv_dt_deepconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%v)
+          status = nf90_inq_varid(ncid, 'dq_dt_deepconv', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%q)
+       endif
+
+       if (iprc == prc_cldMP_cfg(3)) then
+          ! Metadata
+          physics_process(iprc)%order      = iprc
+          physics_process(iprc)%name       = "cldMP"
+          if (prc_cldMP_cfg(1) == 1) then
+             physics_process(iprc)%use_sim = .true.
+          else
+             physics_process(1)%nprg_active = 2
+             iactive_T   = 1
+             iactive_q   = 2
+          endif
+          if (prc_cldMP_cfg(2) == 1) then
+             physics_process(iprc)%time_split = .true.
+          endif
+
+          ! Data
+          status = nf90_inq_varid(ncid, 'dT_dt_micro', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%T)
+          status = nf90_inq_varid(ncid, 'dq_dt_micro', varID)
+          if (status == nf90_noerr) status = nf90_get_var(  ncid, varID, physics_process(iprc)%tend2d%q)
+       endif
+
+       ! Which process-suite is "active"? Is process time-split?
+       if (.not. physics_process(iprc)%use_sim) then
+          physics_process(1)%iactive_scheme = iprc
+          physics_process(1)%active_name    = physics_process(iprc)%name
+          if (physics_process(iprc)%time_split) then
+             physics_process(1)%active_tsp = .true.
+          endif
+       endif
+
+    enddo
+
+    if (physics_process(1)%iactive_scheme == 0) then
+       errflg = 1
+       errmsg = "ERROR: No active suite set for CCPP suite simulator"
+       return
+    endif
+
+    print*, "-----------------------------------"
+    print*, "--- Using CCPP suite simulator ---"
+    print*, "-----------------------------------"
+    do iprc = 1,nprc_sim
+       if (physics_process(iprc)%use_sim) then
+          print*,"  simulate_suite: ", trim(physics_process(iprc)%name)
+          print*,"      order:       ", physics_process(iprc)%order
+          print*,"      time_split:  ", physics_process(iprc)%time_split
+       else
+          print*, "  active_suite:   ", trim(physics_process(1)%active_name)
+          print*, "      order:       ", physics_process(physics_process(1)%iactive_scheme)%order
+          print*, "      time_split : ", physics_process(1)%active_tsp
+       endif
+    enddo
+    print*, "-----------------------------------"
+    print*, "-----------------------------------"
+
+  end subroutine load_ccpp_suite_sim
+
+end module GFS_ccpp_suite_sim_pre
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.meta
new file mode 100644
index 000000000..657d39de0
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_ccpp_suite_sim_pre.meta
@@ -0,0 +1,174 @@
+[ccpp-table-properties]
+  name = GFS_ccpp_suite_sim_pre
+  type = scheme
+  dependencies = ../../hooks/machine.F,module_ccpp_suite_simulator.F90
+
+########################################################################
+[ccpp-arg-table]
+  name = GFS_ccpp_suite_sim_pre_run
+  type = scheme
+[do_ccpp_suite_sim]
+  standard_name = flag_for_ccpp_suite_simulator
+  long_name = flag for ccpp suite simulator
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[physics_process]
+  standard_name = physics_process_type_for_CCPP_suite_simulator
+  long_name = physics process type for CCPP suite simulator
+  units = mixed
+  dimensions = (number_of_physics_process_in_CCPP_suite_simulator)
+  type = base_physics_process
+  intent = in
+[dtend]
+  standard_name = cumulative_change_of_state_variables
+  long_name = diagnostic tendencies for state variables
+  units = mixed
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,cumulative_change_of_state_variables_outer_index_max)
+  type = real
+  kind = kind_phys
+  intent = in
+[dtidx]
+  standard_name = cumulative_change_of_state_variables_outer_index
+  long_name = index of state-variable and process in last dimension of diagnostic tendencies array AKA cumulative_change_index
+  units = index
+  dimensions = (number_of_tracers_plus_one_hundred,number_of_cumulative_change_processes)
+  type = integer
+  intent = in
+[dtp]
+  standard_name = timestep_for_physics
+  long_name = physics timestep
+  units = s
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[index_of_process_dcnv]
+  standard_name = index_of_deep_convection_process_process_in_cumulative_change_index
+  long_name = index of deep convection process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_longwave]
+  standard_name = index_of_longwave_heating_process_in_cumulative_change_index
+  long_name = index of longwave heating process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_shortwave]
+  standard_name = index_of_shortwave_heating_process_in_cumulative_change_index
+  long_name = index of shortwave heating process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_orographic_gwd]
+  standard_name = index_of_orographic_gravity_wave_drag_process_in_cumulative_change_index
+  long_name = index of orographic gravity wave drag process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_scnv]
+  standard_name = index_of_shallow_convection_process_process_in_cumulative_change_index
+  long_name = index of shallow convection process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_pbl]
+  standard_name = index_of_subgrid_scale_vertical_mixing_process_in_cumulative_change_index
+  long_name = index of subgrid scale vertical mixing process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_process_mp]
+  standard_name = index_of_microphysics_process_process_in_cumulative_change_index
+  long_name = index of microphysics transport process in second dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_temperature]
+  standard_name = index_of_temperature_in_cumulative_change_index
+  long_name = index of temperature in first dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_x_wind]
+  standard_name = index_of_x_wind_in_cumulative_change_index
+  long_name = index of x-wind in first dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[index_of_y_wind]
+  standard_name = index_of_y_wind_in_cumulative_change_index
+  long_name = index of x-wind in first dimension of array cumulative change index
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[ntqv]
+  standard_name = index_of_specific_humidity_in_tracer_concentration_array
+  long_name = tracer index for water vapor (specific humidity)
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[active_phys_tend]
+  standard_name = tendencies_for_active_process_in_ccpp_suite_simulator
+  long_name = tendencies for active physics process in ccpp suite simulator
+  units = mixed
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_prognostics_varaibles_in_CCPP_suite_simulator)
+  type = real
+  kind = kind_phys
+  intent = out
+[iactive_T]
+  standard_name = index_for_active_T_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for temperature in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_u]
+  standard_name = index_for_active_u_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for zonal wind in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_v]
+  standard_name = index_for_active_v_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for meridional wind in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_q]
+  standard_name = index_for_active_q_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for moisture in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[errmsg]
+  standard_name = ccpp_error_message
+  long_name = error message for error handling in CCPP
+  units = none
+  dimensions = ()
+  type = character
+  kind = len=*
+  intent = out
+[errflg]
+  standard_name = ccpp_error_code
+  long_name = error code for error handling in CCPP
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = out
\ No newline at end of file
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
index 075bfc039..b338f4ad4 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.F90
@@ -61,7 +61,7 @@ module GFS_phys_time_vary
 !! @{
       subroutine GFS_phys_time_vary_init (                                                         &
               me, master, ntoz, h2o_phys, iaerclm, iccn, iflip, im, nx, ny, idate, xlat_d, xlon_d, &
-              jindx1_o3, jindx2_o3, ddy_o3, ozphys, jindx1_h, jindx2_h, ddy_h, h2opl,fhour,        &
+              jindx1_o3, jindx2_o3, ddy_o3, jindx1_h, jindx2_h, ddy_h, h2opl,fhour,                &
               jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm,            &
               jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap,              &
               do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,                            &
@@ -73,7 +73,7 @@ subroutine GFS_phys_time_vary_init (
               zwtxy, xlaixy, xsaixy, lfmassxy, stmassxy, rtmassxy, woodxy, stblcpxy, fastcpxy,     &
               smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy,   &
               slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, lsm_cold_start, nthrds,    &
-              errmsg, errflg)
+              ozphys, errmsg, errflg)
 
          implicit none
 
@@ -801,7 +801,10 @@ subroutine GFS_phys_time_vary_timestep_init (
                              fhour, iflip, jindx1_aer, jindx2_aer, &
                              ddy_aer, iindx1_aer,           &
                              iindx2_aer, ddx_aer,           &
-                             levs, prsl, aer_nm)
+                             levs, prsl, aer_nm, errmsg, errflg)
+           if(errflg /= 0) then
+             return
+           endif
          endif
          
 !       Not needed for SCM:
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
index a9094a075..e721644b8 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_phys_time_vary.scm.meta
@@ -887,6 +887,13 @@
   dimensions = ()
   type = integer
   intent = in
+[ozphys]
+  standard_name = dataset_for_ozone_physics
+  long_name = dataset for NRL ozone physics
+  units = mixed
+  dimensions = ()
+  type = ty_ozphys
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.F90
new file mode 100644
index 000000000..c1592263d
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.F90
@@ -0,0 +1,212 @@
+! ########################################################################################
+!
+! Description: This suite simulates the evolution of the internal physics state 
+!              represented by a CCPP Suite Definition File (SDF).
+!
+! To activate this suite it must be a) embedded within the SDF and b) activated through
+! the physics namelist.
+! The derived-data type "base_physics_process" contains the metadata needed to reconstruct 
+! the temporal evolution of the state. An array of base_physics_process, physics_process,
+! is populated by the host during initialization and passed to the physics. Additionally, 
+! this type holds any data, or type-bound procedures, required by the suite simulator(s).
+!
+! For this initial demonstration we are using 2-dimensional (height, time) forcing data,
+! which is on the same native vertical grid as the SCM. The dataset has a temporal 
+! resolution of 1-hour, created by averaging all local times from a Tropical Warm Pool 
+! International Cloud Experiment (TWPICE) case. This was to create a dataset with a 
+! (constant) diurnal cycle.
+!
+! ########################################################################################
+module ccpp_suite_simulator
+  use machine, only: kind_phys
+  use module_ccpp_suite_simulator, only: base_physics_process, sim_LWRAD, sim_SWRAD,    &
+       sim_PBL, sim_GWD, sim_DCNV, sim_SCNV, sim_cldMP
+  implicit none
+  public ccpp_suite_simulator_run
+contains
+
+  ! ######################################################################################
+  !
+  ! SUBROUTINE ccpp_suite_simulator_run
+  !
+  ! ######################################################################################
+!! \section arg_table_ccpp_suite_simulator_run
+!! \htmlinclude ccpp_suite_simulator_run.html
+!!
+  subroutine ccpp_suite_simulator_run(do_ccpp_suite_sim, kdt, nCol, nLay, dtp, jdat,     &
+       iactive_T, iactive_u, iactive_v, iactive_q, proc_start, proc_end, physics_process,&
+       in_pre_active, in_post_active, tgrs, ugrs, vgrs, qgrs, active_phys_tend, gt0, gu0,&
+       gv0, gq0, errmsg, errflg)
+
+    ! Inputs
+    logical,           intent(in)  :: do_ccpp_suite_sim
+    integer,           intent(in)  :: kdt, nCol, nLay, jdat(8), iactive_T, iactive_u,    &
+                                      iactive_v, iactive_q
+    real(kind_phys),   intent(in)  :: dtp, tgrs(:,:), ugrs(:,:), vgrs(:,:), qgrs(:,:,:), &
+                                      active_phys_tend(:,:,:)
+    ! Outputs
+    type(base_physics_process),intent(inout) :: physics_process(:)
+    real(kind_phys), intent(inout) :: gt0(:,:), gu0(:,:), gv0(:,:), gq0(:,:)
+    character(len=*),intent(out)   :: errmsg
+    integer,         intent(out)   :: errflg
+    integer,         intent(inout) :: proc_start, proc_end
+    logical,         intent(inout) :: in_pre_active, in_post_active
+
+    ! Locals
+    integer :: iCol, year, month, day, hour, min, sec, iprc
+    real(kind_phys), dimension(nCol,nLay) :: gt1, gu1, gv1, dTdt, dudt, dvdt, gq1, dqdt
+
+    ! Initialize CCPP error handling variables
+    errmsg = ''
+    errflg = 0
+
+    if (.not. do_ccpp_suite_sim) return
+
+    ! Current forecast time (Data-format specific)
+    year  = jdat(1)
+    month = jdat(2)
+    day   = jdat(3)
+    hour  = jdat(5)
+    min   = jdat(6)
+    sec   = jdat(7)
+
+    ! Set state at beginning of the physics timestep.
+    gt1(:,:)  = tgrs(:,:)
+    gu1(:,:)  = ugrs(:,:)
+    gv1(:,:)  = vgrs(:,:)
+    gq1(:,:)  = qgrs(:,:,1)
+    dTdt(:,:) = 0.
+    dudt(:,:) = 0.
+    dvdt(:,:) = 0.
+    dqdt(:,:) = 0.
+
+    !
+    ! Set bookeeping indices
+    !
+    if (in_pre_active) then
+       proc_start = 1
+       proc_end   = max(1,physics_process(1)%iactive_scheme-1)
+    endif
+    if (in_post_active) then
+       proc_start = physics_process(1)%iactive_scheme
+       proc_end   = size(physics_process)
+    endif
+
+    !
+    ! Simulate internal physics timestep evolution.
+    !
+    do iprc = proc_start,proc_end
+       do iCol = 1,nCol
+
+          ! Reset locals
+          physics_process(iprc)%tend1d%T(:) = 0.
+          physics_process(iprc)%tend1d%u(:) = 0.
+          physics_process(iprc)%tend1d%v(:) = 0.
+          physics_process(iprc)%tend1d%q(:) = 0.
+
+          ! Using scheme simulator
+          ! Very simple...
+          ! Interpolate 2D data (time,level) tendency to local time. 
+          ! Here the data is already on the SCM vertical coordinate. 
+          !
+          ! In theory the data can be of any dimensionality and the onus falls on the 
+          ! developer to extend the type "base_physics_process" to work with for their
+          ! application. 
+          !
+          if (physics_process(iprc)%use_sim) then
+             if (physics_process(iprc)%name == "LWRAD") then
+                call sim_LWRAD(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "SWRAD")then
+                call sim_SWRAD(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "GWD")then
+                call sim_GWD(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "PBL")then
+                call sim_PBL(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "SCNV")then
+                call sim_SCNV(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "DCNV")then
+                call sim_DCNV(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+             if (physics_process(iprc)%name == "cldMP")then
+                call sim_cldMP(year, month, day, hour, min, sec, physics_process(iprc))
+             endif
+
+          ! Using data tendency from "active" scheme(s).
+          else
+             if (iactive_T > 0) physics_process(iprc)%tend1d%T = active_phys_tend(iCol,:,iactive_T)
+             if (iactive_u > 0) physics_process(iprc)%tend1d%u = active_phys_tend(iCol,:,iactive_u)
+             if (iactive_v > 0) physics_process(iprc)%tend1d%v = active_phys_tend(iCol,:,iactive_v)
+             if (iactive_q > 0) physics_process(iprc)%tend1d%q = active_phys_tend(iCol,:,iactive_q)
+          endif
+
+          ! Update state now? (time-split scheme)
+          if (physics_process(iprc)%time_split) then
+             gt1(iCol,:)  = gt1(iCol,:) + (dTdt(iCol,:) + physics_process(iprc)%tend1d%T)*dtp
+             gu1(iCol,:)  = gu1(iCol,:) + (dudt(iCol,:) + physics_process(iprc)%tend1d%u)*dtp
+             gv1(iCol,:)  = gv1(iCol,:) + (dvdt(iCol,:) + physics_process(iprc)%tend1d%v)*dtp
+             gq1(iCol,:)  = gq1(iCol,:) + (dqdt(iCol,:) + physics_process(iprc)%tend1d%q)*dtp
+             dTdt(iCol,:) = 0.
+             dudt(iCol,:) = 0.
+             dvdt(iCol,:) = 0.
+             dqdt(iCol,:) = 0.
+          ! Accumulate tendencies, update later? (process-split scheme)
+          else
+             dTdt(iCol,:) = dTdt(iCol,:) + physics_process(iprc)%tend1d%T
+             dudt(iCol,:) = dudt(iCol,:) + physics_process(iprc)%tend1d%u
+             dvdt(iCol,:) = dvdt(iCol,:) + physics_process(iprc)%tend1d%v
+             dqdt(iCol,:) = dqdt(iCol,:) + physics_process(iprc)%tend1d%q
+          endif
+       enddo ! END: Loop over columns
+
+       ! Print diagnostics
+       if (physics_process(iprc)%use_sim) then
+          if (physics_process(iprc)%time_split) then
+             write(*,'(a25,i2,a4,i2,a5,a10,a35)') 'CCPP suite simulator: ',iprc,' of ',proc_end,' ',physics_process(iprc)%name,'time split scheme     (simulated)'
+          else
+             write(*,'(a25,i2,a4,i2,a5,a10,a35)') 'CCPP suite simulator: ',iprc,' of ',proc_end,' ',physics_process(iprc)%name,'process split scheme  (simulated)'
+          endif
+       else
+          if (physics_process(iprc)%time_split) then
+             write(*,'(a25,i2,a4,i2,a5,a10,a35)') 'CCPP suite simulator: ',iprc,' of ',proc_end,' ',physics_process(iprc)%name,'time split scheme        (active)'
+          else
+             write(*,'(a25,i2,a4,i2,a5,a10,a35)') 'CCPP suite simulator: ',iprc,' of ',proc_end,' ',physics_process(iprc)%name,'process split scheme     (active)'
+          endif
+          write(*,'(a25,i2)')                     '       # prog. vars.: ',physics_process(1)%nprg_active
+       endif
+    enddo    ! END: Loop over physics processes
+
+    !
+    ! Update state with accumulated tendencies (process-split only)
+    ! (Suites where active scheme is last physical process)
+    !
+    iprc = minval([iprc,proc_end])
+    if (.not. physics_process(iprc)%time_split) then
+       do iCol = 1,nCol
+          gt0(iCol,:) = gt1(iCol,:) + dTdt(iCol,:)*dtp
+          gu0(iCol,:) = gu1(iCol,:) + dudt(iCol,:)*dtp
+          gv0(iCol,:) = gv1(iCol,:) + dvdt(iCol,:)*dtp
+          gq0(iCol,:) = gq1(iCol,:) + dqdt(iCol,:)*dtp
+       enddo
+    endif
+
+    !
+    ! Update bookeeping indices 
+    !
+    if (in_pre_active) then
+       in_pre_active  = .false.
+       in_post_active = .true.
+    endif
+
+    if (size(physics_process) == proc_end) then
+       in_pre_active  = .true.
+       in_post_active = .false.
+    endif
+
+  end subroutine ccpp_suite_simulator_run
+
+end module ccpp_suite_simulator
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.meta
new file mode 100644
index 000000000..3c91faaeb
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/ccpp_suite_simulator.meta
@@ -0,0 +1,201 @@
+[ccpp-table-properties]
+  name = ccpp_suite_simulator
+  type = scheme
+  dependencies = ../../hooks/machine.F,module_ccpp_suite_simulator.F90
+
+[ccpp-arg-table]
+  name = ccpp_suite_simulator_run
+  type = scheme
+[do_ccpp_suite_sim]
+  standard_name = flag_for_ccpp_suite_simulator
+  long_name = flag for ccpp suite simulator
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[kdt]
+  standard_name = index_of_timestep
+  long_name = current forecast iteration
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[nCol]
+  standard_name = horizontal_loop_extent
+  long_name = horizontal dimension
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[nLay]
+  standard_name = vertical_layer_dimension
+  long_name = number of vertical levels
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[dtp]
+  standard_name = timestep_for_physics
+  long_name = physics timestep
+  units = s
+  dimensions = ()
+  type = real
+  kind = kind_phys
+  intent = in
+[jdat]
+  standard_name = date_and_time_of_forecast_in_united_states_order
+  long_name = current forecast date and time
+  units = none
+  dimensions = (8)
+  type = integer
+  intent = in
+[proc_start]
+  standard_name = index_for_first_physics_process_in_CCPP_suite_simulator
+  long_name = index for first physics process in CCPP suite simulator
+  units     = count
+  dimensions = ()
+  type = integer
+  intent = inout
+[proc_end]
+  standard_name = index_for_last_physics_process_in_CCPP_suite_simulator
+  long_name = index for last physics process in CCPP suite simulator
+  units     = count
+  dimensions = ()
+  type = integer
+  intent = inout
+[in_pre_active]
+  standard_name = flag_to_indicate_location_in_physics_process_loop_before_active_scheme
+  long_name = flag to indicate location in physics process loop before active scheme
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = inout
+[in_post_active]
+  standard_name = flag_to_indicate_location_in_physics_process_loop_after_active_scheme
+  long_name = flag to indicate location in physics process loop after active scheme
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = inout
+[tgrs]
+  standard_name = air_temperature
+  long_name = model layer mean temperature
+  units = K
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[ugrs]
+  standard_name = x_wind
+  long_name = zonal wind
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[vgrs]
+  standard_name = y_wind
+  long_name = meridional wind
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[qgrs]
+  standard_name = tracer_concentration
+  long_name = model layer mean tracer concentration
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_tracers)
+  type = real
+  kind = kind_phys
+  intent = in
+[active_phys_tend]
+  standard_name = tendencies_for_active_process_in_ccpp_suite_simulator
+  long_name = tendencies for active physics process in ccpp suite simulator
+  units = mixed
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,number_of_prognostics_varaibles_in_CCPP_suite_simulator)
+  type = real
+  kind = kind_phys
+  intent = in
+[iactive_T]
+  standard_name = index_for_active_T_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for temperature in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_u]
+  standard_name = index_for_active_u_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for zonal wind in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_v]
+  standard_name = index_for_active_v_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for meridional wind in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[iactive_q]
+  standard_name = index_for_active_q_in_CCPP_suite_simulator
+  long_name = index into active process tracer array for moisture in CCPP suite simulator
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[gt0]
+  standard_name = air_temperature_of_new_state
+  long_name = temperature updated by physics
+  units = K
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
+[gu0]
+  standard_name = x_wind_of_new_state
+  long_name = zonal wind updated by physics
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
+[gv0]
+  standard_name = y_wind_of_new_state
+  long_name = meridional wind updated by physics
+  units = m s-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
+[gq0]
+  standard_name = specific_humidity_of_new_state
+  long_name = tracer concentration updated by physics
+  units = kg kg-1
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = inout
+[physics_process]
+  standard_name = physics_process_type_for_CCPP_suite_simulator
+  long_name = physics process type for CCPP suite simulator
+  units = mixed
+  dimensions = (number_of_physics_process_in_CCPP_suite_simulator)
+  type = base_physics_process
+  intent = inout
+[errmsg]
+  standard_name = ccpp_error_message
+  long_name = error message for error handling in CCPP
+  units = none
+  dimensions = ()
+  type = character
+  kind = len=*
+  intent = out
+[errflg]
+  standard_name = ccpp_error_code
+  long_name = error code for error handling in CCPP
+  units = 1
+  dimensions = ()
+  type = integer
+  intent = out
\ No newline at end of file
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90
new file mode 100644
index 000000000..c4f9fc4e4
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.F90
@@ -0,0 +1,328 @@
+! ########################################################################################
+!
+! This module contains the type, base_physics_process, and supporting subroutines needed
+! by the ccpp suite simulator.
+!
+! ########################################################################################
+module module_ccpp_suite_simulator
+!> \section arg_table_module_ccpp_suite_simulator Argument table
+!! \htmlinclude module_ccpp_suite_simulator.html
+!!
+  use machine, only : kind_phys
+  implicit none
+  
+  public base_physics_process
+
+  ! Type containing 1D (time) physics tendencies. 
+  type phys_tend_1d
+     real(kind_phys), dimension(:), allocatable :: T
+     real(kind_phys), dimension(:), allocatable :: u
+     real(kind_phys), dimension(:), allocatable :: v
+     real(kind_phys), dimension(:), allocatable :: q
+     real(kind_phys), dimension(:), allocatable :: p
+     real(kind_phys), dimension(:), allocatable :: z
+  end type phys_tend_1d
+
+  ! Type containing 2D (lev,time) physics tendencies.
+  type phys_tend_2d
+     real(kind_phys), dimension(:),   allocatable :: time
+     real(kind_phys), dimension(:,:), allocatable :: T
+     real(kind_phys), dimension(:,:), allocatable :: u
+     real(kind_phys), dimension(:,:), allocatable :: v
+     real(kind_phys), dimension(:,:), allocatable :: q
+     real(kind_phys), dimension(:,:), allocatable :: p
+     real(kind_phys), dimension(:,:), allocatable :: z
+  end type phys_tend_2d
+
+  ! Type containing 3D (loc,lev,time) physics tendencies.
+  type phys_tend_3d
+     real(kind_phys), dimension(:),     allocatable :: time
+     real(kind_phys), dimension(:),     allocatable :: lon
+     real(kind_phys), dimension(:),     allocatable :: lat
+     real(kind_phys), dimension(:,:,:), allocatable :: T
+     real(kind_phys), dimension(:,:,:), allocatable :: u
+     real(kind_phys), dimension(:,:,:), allocatable :: v
+     real(kind_phys), dimension(:,:,:), allocatable :: q
+  end type phys_tend_3d
+
+  ! Type containing 4D (lon,lat,lev,time) physics tendencies.
+  type phys_tend_4d
+     real(kind_phys), dimension(:),       allocatable :: time
+     real(kind_phys), dimension(:,:),     allocatable :: lon
+     real(kind_phys), dimension(:,:),     allocatable :: lat
+     real(kind_phys), dimension(:,:,:,:), allocatable :: T
+     real(kind_phys), dimension(:,:,:,:), allocatable :: u
+     real(kind_phys), dimension(:,:,:,:), allocatable :: v
+     real(kind_phys), dimension(:,:,:,:), allocatable :: q
+  end type phys_tend_4d
+
+! This type contains the meta information and data for each physics process.
+
+!> \section arg_table_base_physics_process Argument Table
+!! \htmlinclude base_physics_process.html
+!!
+  type base_physics_process
+     character(len=16)  :: name                 ! Physics process name
+     logical            :: time_split = .false. ! Is process time-split?
+     logical            :: use_sim    = .false. ! Is process "active"?
+     integer            :: order                ! Order of process in process-loop 
+     type(phys_tend_1d) :: tend1d               ! Instantaneous data
+     type(phys_tend_2d) :: tend2d               ! 2-dimensional data
+     type(phys_tend_3d) :: tend3d               ! Not used. Placeholder for 3-dimensional spatial data.
+     type(phys_tend_4d) :: tend4d               ! Not used. Placeholder for 4-dimensional spatio-tempo data.
+     character(len=16)  :: active_name          ! "Active" scheme: Physics process name
+     integer            :: iactive_scheme       ! "Active" scheme: Order of process in process-loop
+     logical            :: active_tsp           ! "Active" scheme: Is process time-split?
+     integer            :: nprg_active          ! "Active" scheme: Number of prognostic variables
+   contains
+     generic,   public  :: linterp => linterp_1D, linterp_2D
+     procedure, private :: linterp_1D
+     procedure, private :: linterp_2D
+     procedure, public  :: find_nearest_loc_2d_1d
+     procedure, public  :: cmp_time_wts
+  end type base_physics_process
+
+contains
+
+  ! ####################################################################################
+  ! Type-bound procedure to compute tendency profile for time-of-day.
+  !
+  ! For use with 1D data (level, time) tendencies with diurnal (24-hr) forcing.
+  ! ####################################################################################
+  function linterp_1D(this, var_name, year, month, day, hour, min, sec) result(err_message)
+    class(base_physics_process), intent(inout) :: this
+    character(len=*), intent(in) :: var_name
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: err_message
+    integer :: ti(1), tf(1), ntime
+    real(kind_phys) :: w1, w2
+
+    ! Interpolation weights
+    call this%cmp_time_wts(year, month, day, hour, min, sec, w1, w2, ti, tf)
+
+    ntime = size(this%tend2d%T(1,:))
+
+    select case(var_name)
+    case("T")
+       if (tf(1) .le. ntime) then
+          this%tend1d%T = w1*this%tend2d%T(:,ti(1)) + w2*this%tend2d%T(:,tf(1))
+       else
+          this%tend1d%T = this%tend2d%T(:,1)
+       endif
+    case("u")
+       if (tf(1) .le. ntime) then
+          this%tend1d%u = w1*this%tend2d%u(:,ti(1)) + w2*this%tend2d%u(:,tf(1))
+       else
+          this%tend1d%u = this%tend2d%u(:,1)
+       endif
+    case("v")
+       if (tf(1) .le. ntime) then
+          this%tend1d%v = w1*this%tend2d%v(:,ti(1)) + w2*this%tend2d%v(:,tf(1))
+       else
+          this%tend1d%v = this%tend2d%v(:,1)
+       endif
+    case("q")
+       if (tf(1) .le. ntime) then
+          this%tend1d%q = w1*this%tend2d%q(:,ti(1)) + w2*this%tend2d%q(:,tf(1))
+       else
+          this%tend1d%q = this%tend2d%q(:,1)
+       endif
+    end select
+
+  end function linterp_1D
+
+  ! ####################################################################################
+  ! Type-bound procedure to compute tendency profile for time-of-day.
+  !
+  ! For use with 2D data (location, level, time) tendencies with diurnal (24-hr) forcing.
+  ! This assumes that the location dimension has a [longitude, latitude] allocated with
+  ! each location.
+  ! ####################################################################################
+  function linterp_2D(this, var_name, lon, lat, year, month, day, hour, min, sec) result(err_message)
+    class(base_physics_process), intent(inout) :: this
+    character(len=*), intent(in) :: var_name
+    integer, intent(in) :: year, month, day, hour, min, sec
+    real(kind_phys), intent(in) :: lon, lat
+    character(len=128) :: err_message
+    integer :: ti(1), tf(1), iNearest
+    real(kind_phys) :: w1, w2
+
+    ! Interpolation weights (temporal)
+    call this%cmp_time_wts(year, month, day, hour, min, sec, w1, w2, ti, tf)
+
+    ! Grab data tendency closest to column [lon,lat]
+    iNearest = this%find_nearest_loc_2d_1d(lon,lat)
+
+    select case(var_name)
+    case("T")
+       this%tend1d%T = w1*this%tend3d%T(iNearest,:,ti(1)) + w2*this%tend3d%T(iNearest,:,tf(1))
+    case("u")
+       this%tend1d%u = w1*this%tend3d%u(iNearest,:,ti(1)) + w2*this%tend3d%u(iNearest,:,tf(1))
+    case("v")
+       this%tend1d%v = w1*this%tend3d%v(iNearest,:,ti(1)) + w2*this%tend3d%v(iNearest,:,tf(1))
+    case("q")
+       this%tend1d%q = w1*this%tend3d%q(iNearest,:,ti(1)) + w2*this%tend3d%q(iNearest,:,tf(1))
+    end select
+  end function linterp_2D
+
+  ! ####################################################################################
+  ! Type-bound procedure to find nearest location.
+  ! For use with linterp_2D, NOT YET IMPLEMENTED.
+  ! ####################################################################################
+  pure function find_nearest_loc_2d_1d(this, lon, lat)
+    class(base_physics_process), intent(in) :: this
+    real(kind_phys), intent(in) :: lon, lat
+    integer :: find_nearest_loc_2d_1d
+
+    find_nearest_loc_2d_1d = 1
+  end function find_nearest_loc_2d_1d
+
+  ! ####################################################################################
+  ! Type-bound procedure to compute linear interpolation weights for a diurnal (24-hour)
+  ! forcing.
+  ! ####################################################################################
+  subroutine cmp_time_wts(this, year, month, day, hour, minute, sec, w1, w2, ti, tf)
+    ! Inputs
+    class(base_physics_process), intent(in) :: this
+    integer, intent(in) :: year, month, day, hour, minute, sec
+    ! Outputs
+    integer,intent(out) :: ti(1), tf(1)
+    real(kind_phys),intent(out) :: w1, w2
+    ! Locals
+    real(kind_phys) :: hrofday
+
+    hrofday = hour*3600. + minute*60. + sec
+    ti = max(hour,1)
+    tf = min(ti + 1,24)
+    w1 = ((hour+1)*3600 - hrofday)/3600
+    w2 = 1 - w1
+
+  end subroutine cmp_time_wts
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_LWRAD( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_LWRAD
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_SWRAD( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_SWRAD
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_GWD( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%u)) then
+       errmsg = process%linterp("u", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%v)) then
+       errmsg = process%linterp("v", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_GWD
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_PBL( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%u)) then
+       errmsg = process%linterp("u", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%v)) then
+       errmsg = process%linterp("v", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%q)) then
+       errmsg = process%linterp("q", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_PBL
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_DCNV( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%u)) then
+       errmsg = process%linterp("u", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%v)) then
+       errmsg = process%linterp("v", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%q)) then
+       errmsg = process%linterp("q", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_DCNV
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_SCNV( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%u)) then
+       errmsg = process%linterp("u", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%v)) then
+       errmsg = process%linterp("v", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%q)) then
+       errmsg = process%linterp("q", year,month,day,hour,min,sec)
+    endif
+
+  end subroutine sim_SCNV
+
+  ! ####################################################################################
+  ! ####################################################################################
+  subroutine sim_cldMP( year, month, day, hour, min, sec, process)
+    type(base_physics_process), intent(inout) :: process
+    integer, intent(in) :: year, month, day, hour, min, sec
+    character(len=128) :: errmsg
+
+    if (allocated(process%tend2d%T)) then
+       errmsg = process%linterp("T", year,month,day,hour,min,sec)
+    endif
+    if (allocated(process%tend2d%q)) then
+       errmsg = process%linterp("q", year,month,day,hour,min,sec)
+    endif
+  end subroutine sim_cldMP
+
+end module module_ccpp_suite_simulator
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.meta
new file mode 100644
index 000000000..55b9e07b1
--- /dev/null
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/module_ccpp_suite_simulator.meta
@@ -0,0 +1,24 @@
+[ccpp-table-properties]
+  name = base_physics_process
+  type = ddt
+  dependencies =
+
+[ccpp-arg-table]
+  name = base_physics_process
+  type = ddt
+
+########################################################################
+[ccpp-table-properties]
+  name = module_ccpp_suite_simulator
+  type = module
+  dependencies = ../../hooks/machine.F
+
+[ccpp-arg-table]
+  name = module_ccpp_suite_simulator
+  type = module
+[base_physics_process]
+  standard_name = base_physics_process
+  long_name = definition of type base_physics_process
+  units = DDT
+  dimensions = ()
+  type = base_physics_process
diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/scm_sfc_flux_spec.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/scm_sfc_flux_spec.F90
index 835b468ff..0196c487c 100644
--- a/physics/Interstitials/UFS_SCM_NEPTUNE/scm_sfc_flux_spec.F90
+++ b/physics/Interstitials/UFS_SCM_NEPTUNE/scm_sfc_flux_spec.F90
@@ -59,8 +59,11 @@ subroutine scm_sfc_flux_spec_run (im, u1, v1, z1, t1, q1, p1, roughness_length,
     
     integer, intent(in)    :: im, lkm
     integer, intent(inout) :: islmsk(:)
+    !WL
+    integer, intent(inout) :: use_lake_model(:)
+    !WL
     logical, intent(in)    :: cplflx, cplice
-    logical, intent(inout) :: dry(:), icy(:), flag_cice(:), wet(:), use_lake_model(:)
+    logical, intent(inout) :: dry(:), icy(:), flag_cice(:), wet(:)!WL, use_lake_model(:)
     real(kind=kind_phys), intent(in)    :: cp, grav, hvap, rd, fvirt, vonKarman, min_seaice, tgice, min_lakeice
     real(kind=kind_phys), intent(in)    :: u1(:), v1(:), z1(:), t1(:), q1(:), p1(:), roughness_length(:), &
       spec_sh_flux(:), spec_lh_flux(:), exner_inverse(:), T_surf(:), oceanfrac(:), lakefrac(:), lakedepth(:)
@@ -212,12 +215,12 @@ subroutine scm_sfc_flux_spec_run (im, u1, v1, z1, t1, q1, p1, roughness_length,
   do i = 1, im
     if ((wet(i) .or. icy(i)) .and. lakefrac(i) > 0.0_kind_phys) then
       if (lkm == 1 .and. lakefrac(i) >= 0.15 .and. lakedepth(i) > 1.0_kind_phys) then
-        use_lake_model(i) = .true.
+        use_lake_model(i) = 1 !WL:.true.
       else
-        use_lake_model(i) = .false.
+        use_lake_model(i) = 0 !WL:.false.
       endif
     else
-      use_lake_model(i) = .false.
+      use_lake_model(i) = 0 !WL.false.
     endif
   enddo
 !