diff --git a/cime/src/drivers/mct/main/cime_comp_mod.F90 b/cime/src/drivers/mct/main/cime_comp_mod.F90
index b07946bb0c3f..f109d1e799ba 100644
--- a/cime/src/drivers/mct/main/cime_comp_mod.F90
+++ b/cime/src/drivers/mct/main/cime_comp_mod.F90
@@ -1906,7 +1906,7 @@ subroutine cime_init()
     !----------------------------------------------------------
 
     areafact_samegrid = .false.
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     if (.not.samegrid_alo) then
        call shr_sys_abort(subname//' ERROR: samegrid_alo is false - Must run with same atm/ocn/lnd grids when configured for scam iop')
     else
diff --git a/cime/src/drivers/moab/main/cime_comp_mod.F90 b/cime/src/drivers/moab/main/cime_comp_mod.F90
index de38faaf189d..8694d4188256 100644
--- a/cime/src/drivers/moab/main/cime_comp_mod.F90
+++ b/cime/src/drivers/moab/main/cime_comp_mod.F90
@@ -1697,7 +1697,7 @@ subroutine cime_init()
     !----------------------------------------------------------
 
     areafact_samegrid = .false.
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     if (.not.samegrid_alo) then
        call shr_sys_abort(subname//' ERROR: samegrid_alo is false - Must run with same atm/ocn/lnd grids when configured for scam iop')
     else
diff --git a/components/cam/bld/build-namelist b/components/cam/bld/build-namelist
index d494b9c2a3d5..277b2cd109af 100755
--- a/components/cam/bld/build-namelist
+++ b/components/cam/bld/build-namelist
@@ -3988,7 +3988,7 @@ if ($cfg->get('scam')) {
 }
 
 # CAM generates IOP file for SCAM
-if ($cfg->get('camiop')) {
+if ($cfg->get('e3smreplay')) {
     add_default($nl, 'inithist', 'val'=>'CAMIOP');
     add_default($nl, 'ndens', 'val'=>'1,1');
     add_default($nl, 'mfilt', 'val'=>'1,10');
diff --git a/components/cam/bld/config_files/definition.xml b/components/cam/bld/config_files/definition.xml
index 78e76ee5d3ed..6c79a9e102f1 100644
--- a/components/cam/bld/config_files/definition.xml
+++ b/components/cam/bld/config_files/definition.xml
@@ -142,9 +142,9 @@ Modifications that allow perturbation growth testing: 0=off, 1=on.
 Configure CAM for single column mode: 0=off, 1=on.  This option only
 supported for the Eulerian dycore.
 </entry>
-<entry id="camiop" valid_values="0,1" value="0">
-Configure CAM to generate an IOP file that can be used to drive SCAM: 0=no, 1=yes.
-This option only supported for the Eulerian dycore.
+<entry id="e3smreplay" valid_values="0,1" value="0">
+Configure E3SM to generate an IOP file that can be used to drive the SCM: 0=no, 1=yes.
+This option only supported for the Spectral Element dycore.
 </entry>
 <entry id="hgrid" value="">
 Horizontal grid specifier.  The recognized values depend on
diff --git a/components/cam/bld/configure b/components/cam/bld/configure
index cb3a741eac18..1ae95c7877f0 100755
--- a/components/cam/bld/configure
+++ b/components/cam/bld/configure
@@ -167,9 +167,9 @@ OPTIONS
 
   Options relevent to SCAM mode:
 
-     -camiop            Configure CAM to generate an IOP file that can be used to drive SCAM.
-                        This switch only works with the Eulerian dycore.
-     -scam              Compiles model in single column mode.  Only works with Eulerian dycore.
+     -e3smreplay        Configure E3SM to generate an IOP file that can be used to drive the SCM.
+                        This switch only works with the Spectral Element dycore.		
+     -scam              Compiles model in single column mode.  Only works with SE dycore.
 
   CAM parallelization:
 
@@ -309,7 +309,7 @@ GetOptions(
     "cam_bld=s"                 => \$opts{'cam_bld'},
     "cam_exe=s"                 => \$opts{'cam_exe'},
     "cam_exedir=s"              => \$opts{'cam_exedir'},
-    "camiop"                    => \$opts{'camiop'},
+    "e3smreplay"                => \$opts{'e3smreplay'},
     "cc=s"                      => \$opts{'cc'},
     "ccsm_seq"                  => \$opts{'ccsm_seq'},
     "cflags=s"                  => \$opts{'cflags'},
@@ -1223,20 +1223,21 @@ if ($print>=2) { print "CAM single column mode (SCAM): $scam$eol"; }
 
 #-----------------------------------------------------------------------------------------------
 # Generate IOP
-if (defined $opts{'camiop'}) {
-    $cfg_ref->set('camiop', 1);
+if (defined $opts{'e3smreplay'}) {
+    $cfg_ref->set('e3smreplay', 1);
 }
-my $camiop = $cfg_ref->get('camiop') ? "ON" : "OFF";
 
-# The only dycore supported in CAMIOP mode is Spectral Element
-if ($camiop eq 'ON' and $dyn_pkg ne 'se') {
+my $e3smreplay = $cfg_ref->get('e3smreplay') ? "ON" : "OFF";
+
+# The only dycore supported in REPLAY mode is Spectral Element
+if ($e3smreplay eq 'ON'  and $dyn_pkg ne 'se') {
     die <<"EOF";
-**  ERROR: CAMIOP mode only works with Spectral Element dycore.
+**  ERROR: REPLAY mode only works with Spectral Element dycore.
 **         Requested dycore is: $dyn_pkg
 EOF
 }
 
-if ($print>=2) { print "Produce IOP file for SCAM: $camiop$eol"; }
+if ($print>=2) { print "Produce IOP file for SCAM: $e3smreplay$eol"; }
 
 #-----------------------------------------------------------------------------------------------
 # Horizontal grid parameters
@@ -1789,7 +1790,7 @@ $cfg_cppdefs .= " -DMAXPATCH_PFT=numpft+1 -DLSMLAT=1 -DLSMLON=1";
 if ($pergro eq "ON") { $cfg_cppdefs .= " -DPERGRO"; }
 
 # Configure CAM to produce IOP files for SCAM
-if ($camiop eq 'ON') { $cfg_cppdefs .= " -DBFB_CAM_SCAM_IOP"; }
+if ($e3smreplay eq 'ON') { $cfg_cppdefs .= " -DE3SM_SCM_REPLAY"; }
 
 # Resolution parameters for rectangular lat/lon grids
 my $nlon = $cfg_ref->get('nlon');
diff --git a/components/cam/src/chemistry/mozart/mo_drydep.F90 b/components/cam/src/chemistry/mozart/mo_drydep.F90
index cc5977a24718..0711fe852a99 100644
--- a/components/cam/src/chemistry/mozart/mo_drydep.F90
+++ b/components/cam/src/chemistry/mozart/mo_drydep.F90
@@ -2036,7 +2036,7 @@ subroutine interp_map( plon, plat, nlon_veg, nlat_veg, npft_veg, lat_veg, lat_ve
                          wetland, vegetation_map, soilw_map, do_soilw )
 
     use mo_constants, only : r2d
-    use scamMod, only : latiop,loniop,scmlat,scmlon,use_camiop
+    use scamMod, only : latiop,loniop,scmlat,scmlon,use_replay
     use shr_scam_mod  , only: shr_scam_getCloseLatLon  ! Standardized system subroutines
     use filenames, only: ncdata
     use dycore, only : dycore_is
@@ -2101,7 +2101,7 @@ subroutine interp_map( plon, plat, nlon_veg, nlat_veg, npft_veg, lat_veg, lat_ve
     ju = plon
 
     if (single_column) then
-       if (use_camiop) then
+       if (use_replay) then
           call getfil (ncdata, ncdata_loc)
           call cam_pio_openfile (piofile, trim(ncdata_loc), PIO_NOWRITE)
           call shr_scam_getCloseLatLon(piofile,scmlat,scmlon,closelat,closelon,latidx,lonidx)
diff --git a/components/cam/src/chemistry/utils/tracer_data.F90 b/components/cam/src/chemistry/utils/tracer_data.F90
index 0d0020baa3c4..fd023e471d46 100644
--- a/components/cam/src/chemistry/utils/tracer_data.F90
+++ b/components/cam/src/chemistry/utils/tracer_data.F90
@@ -1188,7 +1188,7 @@ subroutine read_next_trcdata(flds, file, state )
 
     ! If single column do not interpolate aerosol data, just use the step function. 
     !   The exception is if we are trying to "replay" a column from the full model
-    if(single_column .and. .not. use_camiop) then
+    if(single_column .and. .not. use_replay) then
       file%stepTime = .true.
     endif
 
diff --git a/components/cam/src/control/cam_comp.F90 b/components/cam/src/control/cam_comp.F90
index bd0fe999cdbf..c6424c551734 100644
--- a/components/cam/src/control/cam_comp.F90
+++ b/components/cam/src/control/cam_comp.F90
@@ -87,7 +87,7 @@ subroutine cam_init( cam_out, cam_in, mpicom_atm, &
    use physpkg,          only: phys_init, phys_register
    
    use dycore,           only: dycore_is
-#if (defined BFB_CAM_SCAM_IOP)
+#if (defined E3SM_SCM_REPLAY)
    use history_defaults, only: initialize_iop_history
 #endif
 !   use shr_orb_mod,      only: shr_orb_params
@@ -169,7 +169,7 @@ subroutine cam_init( cam_out, cam_in, mpicom_atm, &
      ! Commented out the hub2atm_alloc call as it overwrite cam_in, which is undesirable. The fields in cam_in are necessary for getting BFB restarts
 	 ! There are no side effects of commenting out this call as this call allocates cam_in and cam_in allocation has already been done in cam_init
      !call hub2atm_alloc( cam_in )
-#if (defined BFB_CAM_SCAM_IOP)
+#if (defined E3SM_SCM_REPLAY)
       call initialize_iop_history()
 #endif
    end if
diff --git a/components/cam/src/control/cam_history.F90 b/components/cam/src/control/cam_history.F90
index 455b604f06d4..ccb60d81ab49 100644
--- a/components/cam/src/control/cam_history.F90
+++ b/components/cam/src/control/cam_history.F90
@@ -3408,7 +3408,7 @@ subroutine h_inquire (t)
       ierr=pio_inq_varid (tape(t)%File,'time_bnds',   tape(t)%tbndid)
       ierr=pio_inq_varid (tape(t)%File,'date_written',tape(t)%date_writtenid)
       ierr=pio_inq_varid (tape(t)%File,'time_written',tape(t)%time_writtenid)
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
       ierr=pio_inq_varid (tape(t)%File,'tsec    ',tape(t)%tsecid)
       ierr=pio_inq_varid (tape(t)%File,'bdate   ',tape(t)%bdateid)
 #endif
@@ -3761,8 +3761,8 @@ subroutine h_define (t, restart)
     str = 'CF-1.0'
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'Conventions', trim(str))
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'source', 'CAM')
-#if ( defined BFB_CAM_SCAM_IOP )
-    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'CAM_GENERATED_FORCING','create SCAM IOP dataset')
+#if defined (E3SM_SCM_REPLAY)
+    ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'E3SM_GENERATED_FORCING','create SCM IOP dataset')
 #endif
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'case',caseid)
     ierr=pio_put_att (tape(t)%File, PIO_GLOBAL, 'title',ctitle)
@@ -3822,7 +3822,7 @@ subroutine h_define (t, restart)
       str = 'base date (YYYYMMDD)'
       ierr=pio_put_att (tape(t)%File, tape(t)%nbdateid, 'long_name', trim(str))
 
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
       ierr=pio_def_var (tape(t)%File,'bdate',PIO_INT,tape(t)%bdateid)
       str = 'base date (YYYYMMDD)'
       ierr=pio_put_att (tape(t)%File, tape(t)%bdateid, 'long_name', trim(str))
@@ -3899,7 +3899,7 @@ subroutine h_define (t, restart)
       end if
 
 
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
       ierr=pio_def_var (tape(t)%File,'tsec ',pio_int,(/timdim/), tape(t)%tsecid)
       str = 'current seconds of current date needed for scam'
       ierr=pio_put_att (tape(t)%File, tape(t)%tsecid, 'long_name', trim(str))
@@ -4149,7 +4149,7 @@ subroutine h_define (t, restart)
 
       ierr = pio_put_var(tape(t)%File, tape(t)%nbdateid, (/nbdate/))
       call cam_pio_handle_error(ierr, 'h_define: cannot put nbdate')
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
       ierr = pio_put_var(tape(t)%File, tape(t)%bdateid, (/nbdate/))
       call cam_pio_handle_error(ierr, 'h_define: cannot put bdate')
 #endif
@@ -4530,7 +4530,7 @@ subroutine wshist (rgnht_in)
     character(len=max_string_len) :: fname ! Filename
     logical :: prev              ! Label file with previous date rather than current
     integer :: ierr
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
     integer :: tsec             ! day component of current time
     integer :: dtime            ! seconds component of current time
 #endif
@@ -4674,7 +4674,7 @@ subroutine wshist (rgnht_in)
           end if
 
           ierr = pio_put_var (tape(t)%File, tape(t)%datesecid,(/start/),(/count1/),(/ncsec/))
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
           dtime = get_step_size()
           tsec=dtime*nstep
           ierr = pio_put_var (tape(t)%File, tape(t)%tsecid,(/start/),(/count1/),(/tsec/))
diff --git a/components/cam/src/control/cam_history_support.F90 b/components/cam/src/control/cam_history_support.F90
index 004cc7df3399..47888edcec38 100644
--- a/components/cam/src/control/cam_history_support.F90
+++ b/components/cam/src/control/cam_history_support.F90
@@ -196,7 +196,7 @@ module cam_history_support
     type(var_desc_t) :: f12vmrid         ! var id for f12 volume mixing ratio
     type(var_desc_t) :: sol_tsiid        ! var id for total solar irradiance (W/m2)
     type(var_desc_t) :: datesecid        ! var id for curent seconds of current date
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
     type(var_desc_t) :: bdateid         ! var id for base date
     type(var_desc_t) :: tsecid        ! var id for curent seconds of current date
 #endif
diff --git a/components/cam/src/control/getinterpnetcdfdata.F90 b/components/cam/src/control/getinterpnetcdfdata.F90
index f9ddd42d16b0..c2dcfb02f7e1 100644
--- a/components/cam/src/control/getinterpnetcdfdata.F90
+++ b/components/cam/src/control/getinterpnetcdfdata.F90
@@ -189,7 +189,7 @@ subroutine getinterpncdata( NCID, camlat, camlon, TimeIdx, &
       outdata(1) = tmp(1)
       return                 ! no need to do interpolation 
    endif
-!   if ( use_camiop .and. nlev.eq.plev) then
+!   if ( use_replay .and. nlev.eq.plev) then
    if ( nlev.eq.plev .or. nlev.eq.plev+1) then
       outData(:nlev)= tmp(:nlev)! no need to do interpolation 
    else
diff --git a/components/cam/src/control/history_defaults.F90 b/components/cam/src/control/history_defaults.F90
index db41aa0f0331..a4a8200aafb4 100644
--- a/components/cam/src/control/history_defaults.F90
+++ b/components/cam/src/control/history_defaults.F90
@@ -22,7 +22,7 @@ module history_defaults
 
   public :: bldfld
 
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
   public :: initialize_iop_history
 #endif
 
@@ -60,7 +60,7 @@ subroutine bldfld ()
   end subroutine bldfld
 
 !#######################################################################
-#if ( defined BFB_CAM_SCAM_IOP  )
+#if ( defined E3SM_SCM_REPLAY  )
   subroutine initialize_iop_history()
 !
 ! !DESCRIPTION: 
diff --git a/components/cam/src/control/scamMod.F90 b/components/cam/src/control/scamMod.F90
index bf5b0ccff5d4..7fd5f5871983 100644
--- a/components/cam/src/control/scamMod.F90
+++ b/components/cam/src/control/scamMod.F90
@@ -136,7 +136,7 @@ module scamMod
                                                ! mo_drydep algorithm
 					       
   real(r8), public ::  scm_relaxation_low      ! lowest level to apply relaxation
-  real(r8), public ::  scm_relaxation_high     ! highest level to apply relaxation					       
+  real(r8), public ::  scm_relaxation_high     ! highest level to apply relaxation
 					       
   real(r8), public, pointer :: loniop(:)
   real(r8), public, pointer :: latiop(:)
@@ -186,7 +186,7 @@ module scamMod
   logical*4, public ::  swrad_off     ! turn off SW radiation (assume night)
   logical*4, public ::  lwrad_off     ! turn off LW radiation
   logical*4, public ::  precip_off    ! turn off precipitation processes
-  logical*4, public ::  use_camiop    ! use cam generated forcing 
+  logical*4, public ::  use_replay    ! use e3sm generated forcing 
   logical*4, public ::  use_3dfrc     ! use 3d forcing
   logical*4, public ::  have_heat_glob ! dataset contains global energy fixer
 
@@ -331,13 +331,13 @@ subroutine scam_setopts( scmlat_in, scmlon_in,iopfile_in,single_column_in, &
         endif
         call wrap_open (iopfile, NF90_NOWRITE, ncid)
 	
-	if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', attnum=i ).EQ. NF90_NOERR ) then
-           use_camiop = .true.
+	if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'E3SM_GENERATED_FORCING', attnum=i ).EQ. NF90_NOERR ) then
+           use_replay = .true.
         else
-           use_camiop = .false.
+           use_replay = .false.
         endif
 	
-	if (dycore_is('SE') .and. use_camiop) then
+	if (dycore_is('SE') .and. use_replay) then
 	  call wrap_inq_dimid( ncid, 'ncol', londimid   )
 	  call wrap_inq_dimlen( ncid, londimid, lonsiz   )
 	  latsiz=lonsiz
@@ -379,7 +379,7 @@ subroutine scam_setopts( scmlat_in, scmlon_in,iopfile_in,single_column_in, &
                  scmlon=ioplon
                  write(iulog,*)'For CAM Generated IOP using closest dataset lat and lon'
               else
-                 if (use_camiop) then
+                 if (use_replay) then
                     call shr_scam_GetCloseLatLon(ncid,scmlat,scmlon,ioplat,ioplon,latidx,lonidx)
                     scmlat=ioplat
                     scmlon=ioplon
@@ -606,7 +606,7 @@ subroutine setiopupdate
 
 end subroutine setiopupdate
 
-  subroutine readiopdata(iop_update_surface,hyam,hybm)
+  subroutine readiopdata(iop_update_phase1,hyam,hybm)
 
 !-----------------------------------------------------------------------
 !     
@@ -631,7 +631,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
 #endif
 !------------------------------Locals-----------------------------------
 !     
-   logical, intent(in) :: iop_update_surface
+   logical, intent(in) :: iop_update_phase1
    integer NCID, status
    integer time_dimID, lev_dimID,lev_varID,mod_dimID,&
            mod_varID,sps_varID,sps_dimID
@@ -678,13 +678,13 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
        'readiopdata.F90', __LINE__)
 
 !
-!     if the dataset is a CAM generated dataset set use_camiop to true
-!       CAM IOP datasets have a global attribute called CAM_GENERATED_IOP      
+!     if the dataset is a CAM generated dataset set use_replay to true
+!       E3SM IOP datasets have a global attribute called E3SM_GENERATED_IOP      
 !
-   if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING',attnum=i ).EQ. NF90_NOERR ) then
-      use_camiop = .true.
+   if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'E3SM_GENERATED_FORCING',attnum=i ).EQ. NF90_NOERR ) then
+      use_replay = .true.
    else
-      use_camiop = .false.
+      use_replay = .false.
    endif
 
 !=====================================================================
@@ -752,7 +752,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
 ! =====================================================
 !     read observed aersol data
 
- if(scm_observed_aero .and. .not. iop_update_surface) then
+ if(scm_observed_aero .and. .not. iop_update_phase1) then
    status = NF90_INQ_DIMID( ncid, 'mod', mod_dimID )
    if ( status .ne. nf90_noerr ) then
       write(iulog,* )'ERROR - readiopdata.F:Could not find variable dim ID  for lev'
@@ -870,7 +870,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
 !
 !CAM generated forcing already has pressure on millibars
 !
-   if (.not. use_camiop) then
+   if (.not. use_replay) then
 !
 !     convert pressure to millibars ( lev is expressed in pascals in iop
 !     datasets )
@@ -892,7 +892,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
    call wrap_inq_dimid(ncid, 'lev', levid)
    call wrap_inq_dimid(ncid, 'time', timeid)
 
-   if (dycore_is('SE') .and. use_camiop) then
+   if (dycore_is('SE') .and. use_replay) then
      strt4(1) = closelonidx
      strt4(2) = iopTimeIdx
      strt4(3) = 1
@@ -908,7 +908,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
    cnt4(3)  = 1
    cnt4(4)  = 1
 
-   if (.not. iop_update_surface) then
+   if (.not. iop_update_phase1) then
 
      status = nf90_inq_varid( ncid, 'Ps', varid   )
      if ( status .ne. nf90_noerr ) then
@@ -931,7 +931,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
 !  the dataset pressure levels to the current
 !  scam model levels
         
-     if ( use_camiop ) then
+     if ( use_replay ) then
        do i = 1, plev
          dplevs( i ) = 1000.0_r8 * hyam( i ) + psobs * hybm( i ) / 100.0_r8
        end do
@@ -950,7 +950,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
          dplevs(i) = psobs/100.0_r8
        end if
      end do
-     if (.not. use_camiop ) then
+     if (.not. use_replay ) then
        nlev = total_levs
      endif
      if ( nlev .eq. 1 ) then
@@ -975,7 +975,7 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
 !
 
 !!!!!!!force fill_end to be .true in getinterpncdata () for temperature !!!!!!!!
-     if ( use_camiop ) then
+     if ( use_replay ) then
        call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'t', have_tsair, &
           tsair(1), .true. , scm_crm_mode, &
           dplevs, nlev, psobs, hyam, hybm, tobs, status )
@@ -999,6 +999,26 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
      else
        have_t = .true.
      endif
+     
+     ! If using REPLAY need to be sure that surface temperature is read in 
+     !  for first radiation call, to ensure b4b (or close) reproducibility. 
+     !  Else, for other SCM cases, it is fine initialize as a cold start.   
+     if (is_first_step() .and. use_replay) then      
+       status = nf90_inq_varid( ncid, 'Tg', varid   )
+       if (status .ne. nf90_noerr) then
+         write(iulog,*)'Could not find variable Tg'
+         if ( have_tsair ) then
+           write(iulog,*) 'UsingTsair'
+           tground = tsair     ! use surface value from T field
+         else
+           write(iulog,*) 'UsingTat lowest level'
+           tground = tobs(plev) 
+         endif
+       else
+         call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tground)
+         have_tg = .true.
+       endif    
+     endif      
 
      status = nf90_inq_varid( ncid, 'qsrf', varid   )
      if ( status .ne. nf90_noerr ) then
@@ -1286,22 +1306,6 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
        ptend= srf(1)
      endif
 
-     call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
-       'omega', .true., ptend, fill_ends, scm_crm_mode, &
-       dplevs, nlev,psobs, hyam, hybm, wfld, status )
-     if ( status .ne. nf90_noerr ) then
-       have_omega = .false.
-       write(iulog,*)'Could not find variable omega'
-       if ( .not. use_userdata ) then
-         status = nf90_close( ncid )
-         return
-       else
-         write(iulog,*) 'Using value from Analysis Dataset'
-       endif
-     else
-       have_omega = .true.
-     endif
-
      call plevs0(1    ,plon   ,plev    ,psobs   ,pint,pmid ,pdel)
      call shr_sys_flush( iulog )
 !
@@ -1456,9 +1460,25 @@ subroutine readiopdata(iop_update_surface,hyam,hybm)
        have_tg = .true.
      endif
      
-     ! If camiop is used, then need to read in the global
+     call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
+       'omega', .true., ptend, fill_ends, scm_crm_mode, &
+        dplevs, nlev,psobs, hyam, hybm, wfld, status )
+     if ( status .ne. nf90_noerr ) then
+        have_omega = .false.
+        write(iulog,*)'Could not find variable omega'
+        if ( .not. use_userdata ) then
+          status = nf90_close( ncid )
+          return
+        else
+          write(iulog,*) 'Using value from Analysis Dataset'
+        endif
+     else
+        have_omega = .true.
+     endif     
+     
+     ! If REPLAY is used, then need to read in the global
      !   energy fixer
-     if (use_camiop) then 
+     if (use_replay) then 
        status = nf90_inq_varid( ncid, 'heat_glob', varid   )
        if (status .ne. nf90_noerr) then
          have_heat_glob = .false.
diff --git a/components/cam/src/control/startup_initialconds.F90 b/components/cam/src/control/startup_initialconds.F90
index fad09cddc685..1fceb90744a7 100644
--- a/components/cam/src/control/startup_initialconds.F90
+++ b/components/cam/src/control/startup_initialconds.F90
@@ -23,7 +23,7 @@ subroutine initial_conds(dyn_in)
 
    use comsrf,        only: initialize_comsrf
    use radae,         only: initialize_radbuffer
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
    use history_defaults, only: initialize_iop_history
 #endif
    
@@ -44,7 +44,7 @@ subroutine initial_conds(dyn_in)
    call initialize_comsrf
    call initialize_radbuffer
 
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
    call initialize_iop_history
 #endif
 
diff --git a/components/cam/src/cpl/atm_comp_mct.F90 b/components/cam/src/cpl/atm_comp_mct.F90
index 0a3dae2e6734..3dc2d347272c 100644
--- a/components/cam/src/cpl/atm_comp_mct.F90
+++ b/components/cam/src/cpl/atm_comp_mct.F90
@@ -51,7 +51,7 @@ module atm_comp_mct
   use cam_logfile      , only: iulog
   use co2_cycle        , only: co2_readFlux_ocn, co2_readFlux_fuel
   use runtime_opts     , only: read_namelist
-  use scamMod          , only: use_camiop,single_column,scmlat,scmlon
+  use scamMod          , only: single_column,scmlat,scmlon
 
 !
 ! !PUBLIC TYPES:
@@ -390,9 +390,6 @@ subroutine atm_init_mct( EClock, cdata_a, x2a_a, a2x_a, NLFilename )
        call seq_timemgr_EClockGetData(EClock,curr_ymd=CurrentYMD, StepNo=StepNo, dtime=DTime_Sync )
        if (StepNo == 0) then
           call atm_import( x2a_a%rattr, cam_in )
-	  if (single_column .and. use_camiop) then
-	    call scam_use_iop_srf( cam_in )
-	  endif
           call cam_run1 ( cam_in, cam_out ) 
           call atm_export( cam_out, a2x_a%rattr )
        else
diff --git a/components/cam/src/dynamics/eul/eul_single_column_mod.F90 b/components/cam/src/dynamics/eul/eul_single_column_mod.F90
index 57d850be3e28..574538f6f709 100644
--- a/components/cam/src/dynamics/eul/eul_single_column_mod.F90
+++ b/components/cam/src/dynamics/eul/eul_single_column_mod.F90
@@ -25,7 +25,7 @@ subroutine scm_setinitial
   integer i, j, k, thelev
   integer inumliq, inumice, icldliq, icldice
 
-  if (.not. use_camiop) then
+  if (.not. use_replay) then
     call cnst_get_ind('NUMLIQ', inumliq, abort=.false.)
     call cnst_get_ind('NUMICE', inumice, abort=.false.)
     call cnst_get_ind('CLDLIQ', icldliq)
diff --git a/components/cam/src/dynamics/eul/inidat.F90 b/components/cam/src/dynamics/eul/inidat.F90
index eca74760573e..fcb7ca5a1315 100644
--- a/components/cam/src/dynamics/eul/inidat.F90
+++ b/components/cam/src/dynamics/eul/inidat.F90
@@ -21,7 +21,7 @@ module inidat
    use spmd_utils,          only: masterproc, mpicom, mpir8
    use cam_control_mod,     only: ideal_phys, aqua_planet, moist_physics, adiabatic
    use cam_initfiles,       only: initial_file_get_id, topo_file_get_id
-   use scamMod,             only: single_column, use_camiop, have_u, have_v, &
+   use scamMod,             only: single_column, use_replay, have_u, have_v, &
                                   have_cldliq, have_cldice,loniop,latiop,scmlat,scmlon
    use cam_logfile,         only: iulog
    use pio,                 only: file_desc_t, pio_noerr, pio_inq_varid, pio_get_att, &
@@ -233,7 +233,7 @@ subroutine read_inidat(fh_ini, fh_topo, dyn_in)
     deallocate ( phis_tmp )
 
     if (single_column) then
-       if ( use_camiop ) then
+       if ( use_replay ) then
           fieldname = 'CLAT1'
           call infld(fieldname, fh_ini, 'lon', 'lat', 1, pcols, begchunk, endchunk, &
                clat2d, readvar, gridname='physgrid')
diff --git a/components/cam/src/dynamics/se/se_single_column_mod.F90 b/components/cam/src/dynamics/se/se_single_column_mod.F90
index dd96da4d888c..d6b1b1d837c3 100644
--- a/components/cam/src/dynamics/se/se_single_column_mod.F90
+++ b/components/cam/src/dynamics/se/se_single_column_mod.F90
@@ -31,7 +31,7 @@ subroutine scm_setinitial(elem)
   integer i, j, k, ie, thelev
   integer inumliq, inumice, icldliq, icldice
 
-  if (.not. use_camiop .and. get_nstep() .eq. 0) then
+  if (.not. use_replay .and. get_nstep() .eq. 0) then
     call cnst_get_ind('NUMLIQ', inumliq, abort=.false.)
     call cnst_get_ind('NUMICE', inumice, abort=.false.)
     call cnst_get_ind('CLDLIQ', icldliq)
@@ -90,10 +90,11 @@ subroutine scm_setinitial(elem)
 
 end subroutine scm_setinitial
 
-subroutine scm_setfield(elem)
+subroutine scm_setfield(elem,iop_update_phase1)
 
   implicit none
 
+  logical, intent(in) :: iop_update_phase1
   type(element_t), intent(inout) :: elem(:)
 
 #ifndef MODEL_THETA_L
@@ -101,9 +102,9 @@ subroutine scm_setfield(elem)
   integer i, j, k, ie
 
   do ie=1,nelemd
-    if (have_ps) elem(ie)%state%ps_v(:,:,:) = psobs 
+    if (have_ps .and. .not. iop_update_phase1) elem(ie)%state%ps_v(:,:,:) = psobs 
     do i=1, PLEV
-      if (have_omega) elem(ie)%derived%omega_p(:,:,i)=wfld(i)  !     set t to tobs at first
+      if (have_omega .and. iop_update_phase1) elem(ie)%derived%omega_p(:,:,i)=wfld(i)  !     set t to tobs at first
     end do
   end do
 
diff --git a/components/cam/src/dynamics/se/stepon.F90 b/components/cam/src/dynamics/se/stepon.F90
index 8cf384092d9c..2842ef038baa 100644
--- a/components/cam/src/dynamics/se/stepon.F90
+++ b/components/cam/src/dynamics/se/stepon.F90
@@ -59,7 +59,7 @@ module stepon
 ! !PRIVATE DATA MEMBERS:
 !
 
-  logical :: iop_update_surface
+  logical :: iop_update_phase1
 
   type (derivative_t)   :: deriv           ! derivative struct
   type (quadrature_t)   :: gv,gp           ! quadratures on velocity and pressure grids
@@ -202,9 +202,9 @@ subroutine stepon_run1( dtime_out, phys_state, phys_tend,               &
   end if
   
   if (single_column) then
-    iop_update_surface = .true. 
-    if (doiopupdate) call readiopdata( iop_update_surface,hyam,hybm )
-    call scm_setfield(elem)       
+    iop_update_phase1 = .true. 
+    if (doiopupdate) call readiopdata( iop_update_phase1,hyam,hybm )
+    call scm_setfield(elem,iop_update_phase1)       
   endif 
   
    call t_barrierf('sync_d_p_coupling', mpicom)
@@ -478,22 +478,23 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
    use cam_logfile, only: iulog
    real(r8), intent(in) :: dtime   ! Time-step
    real(r8) :: ftmp_temp(np,np,nlev,nelemd), ftmp_q(np,np,nlev,pcnst,nelemd)
-   real(r8) :: forcing_temp(npsq,nlev), forcing_q(npsq,nlev,pcnst)
    real(r8) :: out_temp(npsq,nlev), out_q(npsq,nlev), out_u(npsq,nlev), &
                out_v(npsq,nlev), out_psv(npsq)  
    real(r8), parameter :: rad2deg = 180.0 / SHR_CONST_PI
-   real(r8), parameter :: fac = 1000._r8	
-   real(r8) :: term1, term2        
+   real(r8), parameter :: fac = 1000._r8	     
    type(cam_out_t),     intent(inout) :: cam_out(:) ! Output from CAM to surface
    type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
    type (dyn_import_t), intent(inout) :: dyn_in  ! Dynamics import container
    type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
    type (element_t), pointer :: elem(:)
    integer :: rc, i, j, k, p, ie, tl_f
+#if defined (E3SM_SCM_REPLAY)
+   real(r8) :: forcing_temp(npsq,nlev), forcing_q(npsq,nlev,pcnst)
+#endif   
    
    elem => dyn_out%elem
    
-#if (defined BFB_CAM_SCAM_IOP)   
+#if (defined E3SM_SCM_REPLAY)   
 
    tl_f = TimeLevel%n0   ! timelevel which was adjusted by physics
    
@@ -509,11 +510,11 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
      
      ! Update IOP properties e.g. omega, divT, divQ
      
-     iop_update_surface = .false. 
+     iop_update_phase1 = .false. 
      if (doiopupdate) then
        call scm_setinitial(elem)
-       call readiopdata(iop_update_surface,hyam,hybm)
-       call scm_setfield(elem)
+       call readiopdata(iop_update_phase1,hyam,hybm)
+       call scm_setfield(elem,iop_update_phase1)
      endif   
 
    endif   
@@ -524,7 +525,7 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
    call t_stopf  ('dyn_run')
    
    ! Update to get tendency 
-#if (defined BFB_CAM_SCAM_IOP) 
+#if (defined E3SM_SCM_REPLAY) 
 
    tl_f = TimeLevel%n0
    
@@ -533,15 +534,18 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
        do j=1,np
          do i=1,np
 	 
+           ! Note that this calculation will not provide b4b results with 
+	   !  an E3SM because the dynamics tendency is not computed in the exact
+	   !  same way as an E3SM run, introducing error with roundoff 
 	   forcing_temp(i+(j-1)*np,k) = (dyn_in%elem(ie)%state%T(i,j,k,tl_f) - &
-	        ftmp_temp(i,j,k,ie))/dtime - dyn_in%elem(ie)%derived%FT(i,j,k)
+	        ftmp_temp(i,j,k,ie))/dtime - dyn_in%elem(ie)%derived%FT(i,j,k)	
            out_temp(i+(j-1)*np,k) = dyn_in%elem(ie)%state%T(i,j,k,tl_f)
 	   out_u(i+(j-1)*np,k) = dyn_in%elem(ie)%state%v(i,j,1,k,tl_f)
 	   out_v(i+(j-1)*np,k) = dyn_in%elem(ie)%state%v(i,j,2,k,tl_f)
 	   out_q(i+(j-1)*np,k) = dyn_in%elem(ie)%state%Q(i,j,k,1)
 	   out_psv(i+(j-1)*np) = dyn_in%elem(ie)%state%ps_v(i,j,tl_f)
 
-	   do p=1,pcnst		
+	   do p=1,pcnst	 	
 	     forcing_q(i+(j-1)*np,k,p) = (dyn_in%elem(ie)%state%Q(i,j,k,p) - &
 	        ftmp_q(i,j,k,p,ie))/dtime
 	   enddo
@@ -550,12 +554,12 @@ subroutine stepon_run3(dtime, cam_out, phys_state, dyn_in, dyn_out)
        enddo
      enddo
      
-     call outfld('divT3d',forcing_temp,npsq,ie)
      call outfld('Ps',out_psv,npsq,ie)
      call outfld('t',out_temp,npsq,ie)
      call outfld('q',out_q,npsq,ie)
      call outfld('u',out_u,npsq,ie)
      call outfld('v',out_v,npsq,ie)
+     call outfld('divT3d',forcing_temp,npsq,ie)
      do p=1,pcnst
        call outfld(trim(cnst_name(p))//'_dten',forcing_q(:,:,p),npsq,ie)   
      enddo
diff --git a/components/cam/src/physics/cam/cam_diagnostics.F90 b/components/cam/src/physics/cam/cam_diagnostics.F90
index 5c3bec782e82..b55b19a203ee 100644
--- a/components/cam/src/physics/cam/cam_diagnostics.F90
+++ b/components/cam/src/physics/cam/cam_diagnostics.F90
@@ -906,7 +906,7 @@ subroutine diag_phys_writeout(state, psl)
 
 
 
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     call outfld('phis    ',state%phis,    pcols,   lchnk     )
 #endif
 
@@ -1064,7 +1064,7 @@ subroutine diag_phys_writeout(state, psl)
 
     call outfld('OMEGA   ',state%omega,    pcols,   lchnk     )
 
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     call outfld('omega   ',state%omega,    pcols,   lchnk     )
 #endif
 
@@ -1481,7 +1481,7 @@ subroutine diag_conv(state, ztodt, pbuf)
    call outfld('PRECLav ', precl, pcols, lchnk )
    call outfld('PRECCav ', precc, pcols, lchnk )
 
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
    call outfld('Prec   ' , prect, pcols, lchnk )
 #endif
 
@@ -1574,11 +1574,12 @@ subroutine diag_surf (cam_in, cam_out, ps, trefmxav, trefmnav )
     call outfld('RHREFHT',   ftem,      pcols, lchnk)
 
 
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     call outfld('shflx   ',cam_in%shf,   pcols,   lchnk)
     call outfld('lhflx   ',cam_in%lhf,   pcols,   lchnk)
     call outfld('trefht  ',cam_in%tref,  pcols,   lchnk)
     call outfld('Tg', cam_in%ts, pcols, lchnk)
+    call outfld('Tsair',cam_in%ts, pcols, lchnk)
 #endif
 !
 ! Ouput ocn and ice fractions
diff --git a/components/cam/src/physics/cam/check_energy.F90 b/components/cam/src/physics/cam/check_energy.F90
index 74636ff0ccd5..b1003e50b7bb 100644
--- a/components/cam/src/physics/cam/check_energy.F90
+++ b/components/cam/src/physics/cam/check_energy.F90
@@ -814,7 +814,7 @@ subroutine check_energy_fix(state, ptend, nstep, eshflx)
 !------------------------------Arguments--------------------------------
 
     use cam_history, only: outfld
-    use scamMod, only: heat_glob_scm, single_column, use_camiop
+    use scamMod, only: heat_glob_scm, single_column, use_replay
 
     type(physics_state), intent(in   ) :: state
     type(physics_ptend), intent(out)   :: ptend
@@ -838,13 +838,13 @@ subroutine check_energy_fix(state, ptend, nstep, eshflx)
     heat_glob = 0._r8
 #endif
 ! add (-) global mean total energy difference as heating
-    if (single_column .and. use_camiop) then
+    if (single_column .and. use_replay) then
       heat_glob = heat_glob_scm(1)
     endif
     ptend%s(:ncol,:pver) = heat_glob
 !!$    write(iulog,*) "chk_fix: heat", state%lchnk, ncol, heat_glob
 
-#if ( defined BFB_CAM_SCAM_IOP )
+#if ( defined E3SM_SCM_REPLAY )
     if (nstep > 0) then
       heat_out(:ncol) = heat_glob
       call outfld('heat_glob',  heat_out(:ncol), pcols, lchnk)
diff --git a/components/cam/src/physics/cam/physpkg.F90 b/components/cam/src/physics/cam/physpkg.F90
index 2368f3f2bb03..940c0a7860ff 100644
--- a/components/cam/src/physics/cam/physpkg.F90
+++ b/components/cam/src/physics/cam/physpkg.F90
@@ -938,7 +938,7 @@ subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d,  cam_in, cam_out)
     use check_energy,   only: check_energy_gmean
 
     use physics_buffer,         only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate
-#if (defined BFB_CAM_SCAM_IOP )
+#if (defined E3SM_SCM_REPLAY )
     use cam_history,    only: outfld
 #endif
     use comsrf,         only: fsns, fsnt, flns, sgh, sgh30, flnt, landm, fsds
diff --git a/components/cice/src/source/ice_grid.F90 b/components/cice/src/source/ice_grid.F90
index e1289d0974c1..d8dbd6f06dcf 100644
--- a/components/cice/src/source/ice_grid.F90
+++ b/components/cice/src/source/ice_grid.F90
@@ -1191,7 +1191,7 @@ subroutine latlongrid
             tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)
 
             if (single_column) then
-               ULAT  (i,j,iblk) = TLAT(i,j,iblk)+(pi/nj)  
+               ULAT  (i,j,iblk) = TLAT(i,j,iblk)+(pi/ni)  
             else
                if (ny_global == 1) then
                   ULAT  (i,j,iblk) = TLAT(i,j,iblk)
diff --git a/components/clm/src/main/ncdio_pio.F90.in b/components/clm/src/main/ncdio_pio.F90.in
index 238624bb0d62..7727a6f380dc 100644
--- a/components/clm/src/main/ncdio_pio.F90.in
+++ b/components/clm/src/main/ncdio_pio.F90.in
@@ -500,6 +500,12 @@ contains
        ier = pio_inq_dimlen(ncid, dimid, ni)
        if (ier == PIO_NOERR) nj = 1
     end if
+    
+    ier = pio_inq_dimid (ncid, 'topounit', dimid)
+    if (ier == PIO_NOERR) then
+       ier = pio_inq_dimlen(ncid, dimid, ni)
+       if (ier == PIO_NOERR) nj = 1
+    end if       
 
     call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR)
 
@@ -2026,7 +2032,7 @@ end subroutine ncd_io_{DIMS}d_{TYPE}_glob
        if ( trim(dimname)=='nj'.or. trim(dimname)=='lat'.or. trim(dimname)=='lsmlat') then
           start(i)=latidx
           count(i)=1
-       else if ( trim(dimname)=='ni'.or. trim(dimname)=='lon'.or. trim(dimname)=='lsmlon' .or. trim(dimname)=='gridcell') then
+       else if ( trim(dimname)=='ni'.or. trim(dimname)=='lon'.or. trim(dimname)=='lsmlon' .or. trim(dimname)=='gridcell' .or. trim(dimname)=='topounit') then
           start(i)=lonidx
           count(i)=1
        else if ( trim(dimname)=='column') then