diff --git a/Source/Initialization/REMORA_make_new_level.cpp b/Source/Initialization/REMORA_make_new_level.cpp
index 879270f..867189f 100644
--- a/Source/Initialization/REMORA_make_new_level.cpp
+++ b/Source/Initialization/REMORA_make_new_level.cpp
@@ -478,6 +478,14 @@ void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping&
     vec_msku[lev]->setVal(1.0_rt);
     vec_mskv[lev]->setVal(1.0_rt);
 
+    // Initialize these vars even if we aren't using GLS to
+    // avoid issues on e.g. checkpoint
+    vec_tke[lev]->setVal(solverChoice.gls_Kmin);
+    vec_gls[lev]->setVal(solverChoice.gls_Pmin);
+    vec_Lscale[lev]->setVal(0.0_rt);
+    vec_Akk[lev]->setVal(solverChoice.Akk_bak);
+    vec_Akp[lev]->setVal(solverChoice.Akp_bak);
+
     // NOTE: Used to set vec_pm and vec_pn to 1e34 here to make foextrap work
     // when init_type = real. However, this does not appear to be necessary so removing
 
diff --git a/Source/REMORA.H b/Source/REMORA.H
index 24b6cd9..5da2e0d 100644
--- a/Source/REMORA.H
+++ b/Source/REMORA.H
@@ -607,6 +607,7 @@ public:
     void stretch_transform (int lev);
 
     void init_set_vmix (int lev);
+    void set_analytical_vmix (int lev);
     void init_gls_vmix (int lev, SolverChoice solver_choice);
     void gls_prestep (int lev, amrex::MultiFab* mf_gls, amrex::MultiFab* mf_tke,
                       amrex::MultiFab& mf_W,
diff --git a/Source/REMORA.cpp b/Source/REMORA.cpp
index e48c170..bd813e6 100644
--- a/Source/REMORA.cpp
+++ b/Source/REMORA.cpp
@@ -552,11 +552,7 @@ void
 REMORA::init_set_vmix(int lev) {
     Real time = 0.0_rt;
     if (solverChoice.vert_mixing_type == VertMixingType::analytical) {
-        init_custom_vmix(geom[lev], *vec_Akv[lev], *vec_Akt[lev], *vec_z_w[lev], solverChoice);
-        FillPatch(lev, time, *vec_Akv[lev], GetVecOfPtrs(vec_Akv),BdyVars::null,0,true,false);
-        for (int n=0; n<NCONS;n++) {
-            FillPatch(lev, time, *vec_Akt[lev], GetVecOfPtrs(vec_Akt),BdyVars::null,0,false,false);
-        }
+        set_analytical_vmix(lev);
     } else if (solverChoice.vert_mixing_type == VertMixingType::GLS) {
         init_gls_vmix(lev, solverChoice);
         // The GLS initialization just sets the multifab to a value, so there's
@@ -566,6 +562,16 @@ REMORA::init_set_vmix(int lev) {
     }
 }
 
+void
+REMORA::set_analytical_vmix(int lev) {
+    Real time = 0.0_rt;
+    init_custom_vmix(geom[lev], *vec_Akv[lev], *vec_Akt[lev], *vec_z_w[lev], solverChoice);
+    FillPatch(lev, time, *vec_Akv[lev], GetVecOfPtrs(vec_Akv),BdyVars::null,0,true,false);
+    for (int n=0; n<NCONS;n++) {
+        FillPatch(lev, time, *vec_Akt[lev], GetVecOfPtrs(vec_Akt),BdyVars::null,0,false,false);
+    }
+}
+
 void
 REMORA::set_hmixcoef(int lev)
 {
diff --git a/Source/TimeIntegration/REMORA_setup_step.cpp b/Source/TimeIntegration/REMORA_setup_step.cpp
index 8395407..26175c9 100644
--- a/Source/TimeIntegration/REMORA_setup_step.cpp
+++ b/Source/TimeIntegration/REMORA_setup_step.cpp
@@ -176,7 +176,7 @@ REMORA::setup_step (int lev, Real time, Real dt_lev)
 
     if (solverChoice.vert_mixing_type == VertMixingType::analytical) {
         // Update Akv if using analytical mixing
-        init_set_vmix(lev);
+        set_analytical_vmix(lev);
     }
 
     set_zeta_to_Ztavg(lev);