diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml
index b42f68d6f7e9..e522e7412c7b 100644
--- a/components/mpas-albany-landice/src/Registry.xml
+++ b/components/mpas-albany-landice/src/Registry.xml
@@ -492,6 +492,10 @@
                             description="Order of Runge-Kutta time integration to use. A value of 1 is equivalent to forward euler. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods."
                             possible_values="1, 2, 3"
                 />
+		<nml_option name="config_rk3_stages" type="integer" default_value="3" units="stages"
+                        description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2."
+                        possible_values="3 or 4"
+		/>
 		<nml_option name="config_adaptive_timestep" type="logical" default_value=".false." units="unitless"
 			description="Determines if the time step should be adjusted based on the CFL condition or should be steady in time. If true, the config_dt_* options are ignored."
 			possible_values=".true. or .false."
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
index 527032420369..17714dfb6ffc 100644
--- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
+++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
@@ -112,7 +112,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
       character (len=StrKIND), pointer :: config_thickness_advection
       character (len=StrKIND), pointer :: config_thermal_solver
       character (len=StrKIND), pointer :: config_time_integration
-      integer, pointer :: config_rk_order
+      integer, pointer :: config_rk_order, config_rk3_stages
       integer :: rkStage, iCell, iTracer, k
       real (kind=RKIND), dimension(:,:), pointer :: layerThickness, &
                                                     temperature, &
@@ -133,7 +133,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
       real (kind=RKIND), pointer :: deltat
       real (kind=RKIND) :: deltatFull
       real (kind=RKIND), dimension(4) :: rkSubstepWeights
-      real (kind=RKIND), dimension(3) :: rkSSPweights
+      real (kind=RKIND), dimension(4) :: rkSSPweights
       integer :: nRKstages
 
       err = 0
@@ -145,6 +145,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
       call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection)
       call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver)
       call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order)
+      call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages)
       call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration)
 
       call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
@@ -262,16 +263,33 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
 
          rkSSPweights(1) = 1.0_RKIND
          rkSSPweights(2) = 0.5_RKIND
+         rkSSPweights(3) = 0.0_RKIND
+         rkSSPweights(4) = 0.0_RKIND
       elseif ( (trim(config_time_integration) == 'runge_kutta') .and. &
                (config_rk_order == 3) ) then
-         ! use Strong Stability Preserving RK3
-         nRKstages = 3
-         rkSubstepWeights(:) = 1.0_RKIND
-
-         rkSSPweights(1) = 1.0_RKIND
-         rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND
-         rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND
-        ! Add option for 4-stage SSP-RK3? Allows for CFL=2.
+         if (config_rk3_stages == 3) then
+            ! use Strong Stability Preserving RK3
+            nRKstages = 3
+            rkSubstepWeights(:) = 1.0_RKIND
+
+            rkSSPweights(1) = 1.0_RKIND
+            rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND
+            rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND
+            rkSSPweights(4) = 0.0_RKIND
+        ! 4-stage SSP-RK3? Allows for CFL=2.
+        elseif (config_rk3_stages == 4) then
+            nRKstages = 4
+            rkSubstepWeights(:) = 0.5_RKIND
+
+            rkSSPweights(1) = 1.0_RKIND
+            rkSSPweights(2) = 0.0_RKIND
+            rkSSPweights(3) = 2.0_RKIND / 3.0_RKIND
+            rkSSPweights(4) = 0.0_RKIND
+        else
+            err = 1
+            call mpas_log_write('config_rk3_stages must 3 or 4', &
+                                messageType=MPAS_LOG_ERR)
+        endif
       else
          err = 1
          call mpas_log_write('config_time_integration = ' // trim(config_time_integration) &
@@ -299,8 +317,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
          call mpas_timer_stop("advect thickness and tracers")
 
          ! If using SSP RK, then update thickness and tracers incrementally.
-         ! For first RK stage, thickness and tracer updates above are sufficient
-         if ( rkStage > 1 )  then
+         ! For first RK stage, thickness and tracer updates above are sufficient.
+         ! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update.
+         if ( (rkStage > 1) .and. (rkStage < 4) )  then
             call mpas_pool_get_array(geometryPool, 'thickness', thickness)
             call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
             call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d)
@@ -507,7 +526,7 @@ subroutine prepare_advection(domain, err)
       real (kind=RKIND), dimension(:,:), pointer :: normalVelocity
       real (kind=RKIND), dimension(:,:), pointer :: layerNormalVelocity
       real (kind=RKIND), pointer :: calvingCFLdt, faceMeltingCFLdt
-      integer, pointer :: processLimitingTimestep
+      integer, pointer :: processLimitingTimestep, config_rk_order, config_rk3_stages
       integer, dimension(:), pointer :: edgeMask
 
       logical, pointer :: config_print_thickness_advection_info
@@ -515,7 +534,8 @@ subroutine prepare_advection(domain, err)
       logical, pointer :: config_adaptive_timestep_include_DCFL
 
       character (len=StrKIND), pointer :: &
-           config_thickness_advection   ! method for advecting thickness and tracers
+           config_thickness_advection, &  ! method for advecting thickness and tracers
+           config_time_integration
 
       integer :: &
            allowableAdvecDtProcNumberHere, &
@@ -567,6 +587,9 @@ subroutine prepare_advection(domain, err)
       call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep)
       call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep_include_DCFL', config_adaptive_timestep_include_DCFL)
       call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection)
+      call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration)
+      call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order)
+      call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages)
 
       if (trim(config_thickness_advection) == 'none') then
          if (config_adaptive_timestep) then
@@ -612,7 +635,6 @@ subroutine prepare_advection(domain, err)
          err = ior(err, err_tmp)
 
          allowableAdvecDtOnProc = min(allowableAdvecDtOnProc, allowableAdvecDt)
-
          ! Calculate diffusive CFL timestep, if needed
          ! This used to be only calculated if (config_adaptive_timestep_include_DCFL) but for simplicity,
          ! now it is always calculated.  That allows assessment of the DCFL even when it is not being obeyed
@@ -637,6 +659,12 @@ subroutine prepare_advection(domain, err)
          block => block % next
       end do
 
+      ! If using 4-stage SSPRK3, CFL number of 2 is theoretically allowed
+      if ( (trim(config_time_integration) == 'runge_kutta') .and. &
+         (config_rk_order == 3) .and. (config_rk3_stages == 4) ) then
+          allowableAdvecDtOnProc = allowableAdvecDtOnProc * 2.0_RKIND
+          allowableDiffDtOnProc = allowableDiffDtOnProc * 2.0_RKIND
+      endif
 
       ! Local advective CFL info
       call mpas_set_timeInterval(allowableAdvecDtOnProcInterval, dt=allowableAdvecDtOnProc, ierr=err_tmp)