diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90
index 3a3ba6920c..670d1a2fca 100644
--- a/src/core/MOM_CoriolisAdv.F90
+++ b/src/core/MOM_CoriolisAdv.F90
@@ -25,6 +25,7 @@ module MOM_CoriolisAdv
 
 !> Control structure for mom_coriolisadv
 type, public :: CoriolisAdv_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   integer :: Coriolis_Scheme !< Selects the discretization for the Coriolis terms.
                              !! Valid values are:
                              !! - SADOURNY75_ENERGY - Sadourny, 1975
@@ -247,6 +248,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "MOM_CoriolisAdv: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_CoriolisAdv: Module must be initialized before it is used.")
+
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
   vol_neglect = GV%H_subroundoff * (1e-4 * US%m_to_L)**2
@@ -1131,6 +1136,8 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   CS%diag => diag ; CS%Time => Time
 
   ! Read all relevant parameters and write them to the model log.
diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 3100699e6f..734be4c758 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -34,6 +34,7 @@ module MOM_PressureForce_FV
 
 !> Finite volume pressure gradient control structure
 type, public :: PressureForce_FV_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: tides          !< If true, apply tidal momentum forcing.
   real    :: Rho0           !< The density used in the Boussinesq
                             !! approximation [R ~> kg m-3].
@@ -165,6 +166,10 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
+
   if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//&
        "implemented in non-Boussinesq mode.")
@@ -502,6 +507,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   if (.not.associated(CS)) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")
+
   use_p_atm = associated(p_atm)
   use_EOS = associated(tv%eqn_of_state)
   do i=Isq,Ieq+1 ; p0(i) = 0.0 ; enddo
@@ -820,6 +828,8 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag ; CS%Time => Time
   if (associated(tides_CSp)) CS%tides_CSp => tides_CSp
 
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index 27aaf49276..1aa4d9c384 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -30,6 +30,7 @@ module MOM_PressureForce_Mont
 
 !> Control structure for the Montgomery potential form of pressure gradient
 type, public :: PressureForce_Mont_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: tides          !< If true, apply tidal momentum forcing.
   real    :: Rho0           !< The density used in the Boussinesq
                             !! approximation [R ~> kg m-3].
@@ -139,6 +140,10 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
       "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+      "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
   if (use_EOS) then
     if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, &
       "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
@@ -426,6 +431,10 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
        "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
   if (use_EOS) then
     if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, &
       "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
@@ -839,6 +848,8 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag ; CS%Time => Time
   if (associated(tides_CSp)) CS%tides_CSp => tides_CSp
 
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 219d22cc93..a90deeaf58 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -697,6 +697,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
       "btstep: Module MOM_barotropic must be initialized before it is used.")
+
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "btstep: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -2769,6 +2773,10 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
       "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
+
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   MS%isdw = G%isd ; MS%iedw = G%ied ; MS%jsdw = G%jsd ; MS%jedw = G%jed
@@ -3306,6 +3314,10 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
 ! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
   if (.not.associated(CS)) call MOM_error(FATAL, &
       "btcalc: Module MOM_barotropic must be initialized before it is used.")
+
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "btcalc: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
 
   use_default = .false.
@@ -4198,6 +4210,10 @@ subroutine bt_mass_source(h, eta, set_cor, G, GV, CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "bt_mass_source: "// &
         "Module MOM_barotropic must be initialized before it is used.")
+
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, "bt_mass_source: "// &
+        "Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90
index a9cd01a6df..35965402bb 100644
--- a/src/core/MOM_continuity_PPM.F90
+++ b/src/core/MOM_continuity_PPM.F90
@@ -26,6 +26,7 @@ module MOM_continuity_PPM
 
 !> Control structure for mom_continuity_ppm
 type, public :: continuity_PPM_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   type(diag_ctrl), pointer :: diag !< Diagnostics control structure.
   logical :: upwind_1st      !< If true, use a first-order upwind scheme.
   logical :: monotonic       !< If true, use the Colella & Woodward monotonic
@@ -136,6 +137,10 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "MOM_continuity_PPM: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_continuity_PPM: Module must be initialized before it is used.")
+
   x_first = (MOD(G%first_direction,2) == 0)
 
   if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, &
@@ -2210,6 +2215,8 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS)
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
 ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
   call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", CS%monotonic, &
diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90
index 7817fc4959..428b6f79ae 100644
--- a/src/diagnostics/MOM_diagnostics.F90
+++ b/src/diagnostics/MOM_diagnostics.F90
@@ -49,6 +49,7 @@ module MOM_diagnostics
 
 !> The control structure for the MOM_diagnostics module
 type, public :: diagnostics_CS ; private
+  logical :: initialized = .false.     !< True if this control structure has been initialized.
   real :: mono_N2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
                                        !! monotonic for the purposes of calculating the equivalent
                                        !! barotropic wave speed.
@@ -272,6 +273,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
   if (loc(CS)==0) call MOM_error(FATAL, &
          "calculate_diagnostic_fields: Module must be initialized before used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, &
+         "calculate_diagnostic_fields: Module must be initialized before used.")
+
   call calculate_derivs(dt, G, CS)
 
   if (dt > 0.0) then
@@ -1253,6 +1257,9 @@ subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "register_time_deriv: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "register_time_deriv: Module must be initialized before it is used.")
+
   if (CS%num_time_deriv >= MAX_FIELDS_) then
     call MOM_error(WARNING,"MOM_diagnostics:  Attempted to register more than " // &
                    "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
@@ -1620,6 +1627,8 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   CS%diag => diag
   use_temperature = associated(tv%T)
   call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90
index 602041372b..4c70245a1f 100644
--- a/src/diagnostics/MOM_sum_output.F90
+++ b/src/diagnostics/MOM_sum_output.F90
@@ -59,6 +59,8 @@ module MOM_sum_output
 
 !> The control structure for the MOM_sum_output module
 type, public :: sum_output_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
+
   type(Depth_List)              :: DL !< The sorted depth list.
 
   integer, allocatable, dimension(:) :: lH
@@ -160,6 +162,8 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, &
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
   call get_param(param_file, mdl, "CALCULATE_APE", CS%do_APE_calc, &
@@ -490,6 +494,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "write_energy: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "write_energy: Module must be initialized before it is used.")
+
   do j=js,je ; do i=is,ie
     areaTm(i,j) = G%mask2dT(i,j)*G%areaT(i,j)
   enddo ; enddo
diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90
index 833e7d8165..dd6777d032 100644
--- a/src/diagnostics/MOM_wave_speed.F90
+++ b/src/diagnostics/MOM_wave_speed.F90
@@ -26,6 +26,7 @@ module MOM_wave_speed
 
 !> Control structure for MOM_wave_speed
 type, public :: wave_speed_CS ; private
+  logical :: initialized = .false.     !< True if this control structure has been initialized.
   logical :: use_ebt_mode = .false.    !< If true, calculate the equivalent barotropic wave speed instead
                                        !! of the first baroclinic wave speed.
                                        !! This parameter controls the default behavior of wave_speed() which
@@ -149,6 +150,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_wave_speed: "// &
            "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// &
+           "Module must be initialized before it is used.")
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -731,6 +736,11 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
            "Module must be initialized before it is used.")
   endif
 
+  if (present(CS)) then
+    if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// &
+           "Module must be initialized before it is used.")
+  endif
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -1200,6 +1210,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   ! Write all relevant parameters to the model log.
   call log_version(mdl, version)
 
diff --git a/src/diagnostics/MOM_wave_structure.F90 b/src/diagnostics/MOM_wave_structure.F90
index cf4c518889..41464845ae 100644
--- a/src/diagnostics/MOM_wave_structure.F90
+++ b/src/diagnostics/MOM_wave_structure.F90
@@ -36,6 +36,7 @@ module MOM_wave_structure
 
 !> The control structure for the MOM_wave_structure module
 type, public :: wave_structure_CS ; !private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                    !! regulate the timing of diagnostic output.
   real, allocatable, dimension(:,:,:) :: w_strct
@@ -198,6 +199,9 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
            "Module must be initialized before it is used.")
   !endif
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_structure: "// &
+         "Module must be initialized before it is used.")
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -742,6 +746,8 @@ subroutine wave_structure_init(Time, G, GV, param_file, diag, CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", CS%int_tide_source_x, &
                  "X Location of generation site for internal tide", default=1.)
   call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", CS%int_tide_source_y, &
diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90
index 5d81db10a3..8dcae65d52 100644
--- a/src/framework/MOM_restart.F90
+++ b/src/framework/MOM_restart.F90
@@ -69,6 +69,7 @@ module MOM_restart
 
 !> A restart registry and the control structure for restarts
 type, public :: MOM_restart_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: restart    !< restart is set to .true. if the run has been started from a full restart
                         !! file.  Otherwise some fields must be initialized approximately.
   integer :: novars = 0 !< The number of restart fields that have been registered.
@@ -156,6 +157,9 @@ subroutine register_restart_field_ptr3d(f_ptr, var_desc, mandatory, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "register_restart_field: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "register_restart_field: Module must be initialized before it is used.")
+
   call lock_check(CS, var_desc)
 
   CS%novars = CS%novars+1
@@ -189,6 +193,9 @@ subroutine register_restart_field_ptr4d(f_ptr, var_desc, mandatory, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "register_restart_field: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "register_restart_field: Module must be initialized before it is used.")
+
   call lock_check(CS, var_desc)
 
   CS%novars = CS%novars+1
@@ -222,6 +229,9 @@ subroutine register_restart_field_ptr2d(f_ptr, var_desc, mandatory, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "register_restart_field: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "register_restart_field: Module must be initialized before it is used.")
+
   call lock_check(CS, var_desc)
 
   CS%novars = CS%novars+1
@@ -254,6 +264,9 @@ subroutine register_restart_field_ptr1d(f_ptr, var_desc, mandatory, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "register_restart_field: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "register_restart_field: Module must be initialized before it is used.")
+
   call lock_check(CS, var_desc)
 
   CS%novars = CS%novars+1
@@ -286,6 +299,9 @@ subroutine register_restart_field_ptr0d(f_ptr, var_desc, mandatory, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "register_restart_field: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "register_restart_field: Module must be initialized before it is used.")
+
   call lock_check(CS, var_desc)
 
   CS%novars = CS%novars+1
@@ -397,6 +413,10 @@ subroutine register_restart_field_4d(f_ptr, name, mandatory, CS, longname, units
       "register_restart_field_4d: Module must be initialized before "//&
       "it is used to register "//trim(name))
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // &
+      "register_restart_field_4d: Module must be initialized before "//&
+      "it is used to register "//trim(name))
+
   call lock_check(CS, name=name)
 
   vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, &
@@ -427,6 +447,10 @@ subroutine register_restart_field_3d(f_ptr, name, mandatory, CS, longname, units
       "register_restart_field_3d: Module must be initialized before "//&
       "it is used to register "//trim(name))
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // &
+      "register_restart_field_3d: Module must be initialized before "//&
+      "it is used to register "//trim(name))
+
   call lock_check(CS, name=name)
 
   vd = var_desc(name, units=units, longname=longname, hor_grid=hor_grid, &
@@ -457,6 +481,11 @@ subroutine register_restart_field_2d(f_ptr, name, mandatory, CS, longname, units
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart: " // &
       "register_restart_field_2d: Module must be initialized before "//&
       "it is used to register "//trim(name))
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // &
+      "register_restart_field_2d: Module must be initialized before "//&
+      "it is used to register "//trim(name))
+
   zgrid = '1' ; if (present(z_grid)) zgrid = z_grid
 
   call lock_check(CS, name=name)
@@ -488,6 +517,11 @@ subroutine register_restart_field_1d(f_ptr, name, mandatory, CS, longname, units
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart: " // &
       "register_restart_field_3d: Module must be initialized before "//&
       "it is used to register "//trim(name))
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // &
+      "register_restart_field_3d: Module must be initialized before "//&
+      "it is used to register "//trim(name))
+
   hgrid = '1' ; if (present(hor_grid)) hgrid = hor_grid
 
   call lock_check(CS, name=name)
@@ -517,6 +551,10 @@ subroutine register_restart_field_0d(f_ptr, name, mandatory, CS, longname, units
       "register_restart_field_0d: Module must be initialized before "//&
       "it is used to register "//trim(name))
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart: " // &
+      "register_restart_field_0d: Module must be initialized before "//&
+      "it is used to register "//trim(name))
+
   call lock_check(CS, name=name)
 
   vd = var_desc(name, units=units, longname=longname, hor_grid='1', &
@@ -538,6 +576,10 @@ function query_initialized_name(name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -570,6 +612,10 @@ function query_initialized_0d(f_ptr, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -595,6 +641,10 @@ function query_initialized_1d(f_ptr, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -621,6 +671,10 @@ function query_initialized_2d(f_ptr, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -647,6 +701,10 @@ function query_initialized_3d(f_ptr, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -673,6 +731,10 @@ function query_initialized_4d(f_ptr, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -700,6 +762,10 @@ function query_initialized_0d_name(f_ptr, name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -734,6 +800,10 @@ function query_initialized_1d_name(f_ptr, name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -768,6 +838,10 @@ function query_initialized_2d_name(f_ptr, name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -802,6 +876,10 @@ function query_initialized_3d_name(f_ptr, name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -836,6 +914,10 @@ function query_initialized_4d_name(f_ptr, name, CS) result(query_initialized)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "query_initialized: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "query_initialized: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   query_initialized = .false.
@@ -908,6 +990,10 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "save_restart: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "save_restart: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   ! With parallel read & write, it is possible to disable the following...
@@ -1099,6 +1185,10 @@ subroutine restore_state(filename, directory, day, G, CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "restore_state: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "restore_state: Module must be initialized before it is used.")
+
   if (CS%novars > CS%max_fields) call restart_error(CS)
 
   ! Get NetCDF ids for all of the restart files.
@@ -1295,6 +1385,9 @@ function restart_files_exist(filename, directory, G, CS)
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "restart_files_exist: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "restart_files_exist: Module must be initialized before it is used.")
+
   if ((LEN_TRIM(filename) == 1) .and. (filename(1:1) == 'F')) then
     num_files = get_num_restart_files('r', directory, G, CS)
   else
@@ -1320,6 +1413,9 @@ function determine_is_new_run(filename, directory, G, CS) result(is_new_run)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "determine_is_new_run: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "determine_is_new_run: Module must be initialized before it is used.")
   if (LEN_TRIM(filename) > 1) then
     CS%new_run = .false.
   elseif (LEN_TRIM(filename) == 0) then
@@ -1346,6 +1442,10 @@ function is_new_run(CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "is_new_run: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "is_new_run: Module must be initialized before it is used.")
+
   if (.not.CS%new_run_set) call MOM_error(FATAL, "MOM_restart " // &
       "determine_is_new_run must be called for a restart file before is_new_run.")
 
@@ -1391,6 +1491,9 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths,
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "open_restart_units: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "open_restart_units: Module must be initialized before it is used.")
+
   ! Get NetCDF ids for all of the restart files.
   num_restart = 0 ; nf = 0 ; start_char = 1
   do while (start_char <= len_trim(filename) )
@@ -1506,6 +1609,9 @@ function get_num_restart_files(filenames, directory, G, CS, file_paths) result(n
   if (.not.associated(CS)) call MOM_error(FATAL, "MOM_restart " // &
       "get_num_restart_files: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, "MOM_restart " // &
+      "get_num_restart_files: Module must be initialized before it is used.")
+
   ! This call uses open_restart_units without the optional arguments needed to actually
   ! open the files to determine the number of restart files.
   num_files = open_restart_units(filenames, directory, G, CS, file_paths=file_paths)
@@ -1535,6 +1641,8 @@ subroutine restart_init(param_file, CS, restart_root)
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   ! Determine whether all paramters are set to their default values.
   call get_param(param_file, mdl, "PARALLEL_RESTARTFILES", CS%parallel_restartfiles, &
                  default=.false., do_not_log=.true.)
diff --git a/src/framework/MOM_write_cputime.F90 b/src/framework/MOM_write_cputime.F90
index c9200cf41c..9df994448b 100644
--- a/src/framework/MOM_write_cputime.F90
+++ b/src/framework/MOM_write_cputime.F90
@@ -20,6 +20,7 @@ module MOM_write_cputime
 
 !> A control structure that regulates the writing of CPU time
 type, public :: write_cputime_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real :: maxcpu                !<   The maximum amount of cpu time per processor
                                 !! for which MOM should run before saving a restart
                                 !! file and quiting with a return value that
@@ -71,6 +72,8 @@ subroutine MOM_write_cputime_init(param_file, directory, Input_start_time, CS)
     CS%prev_cputime = new_cputime
   endif
 
+  CS%initialized = .true.
+
   ! Read all relevant parameters and write them to the model log.
 
   ! Determine whether all paramters are set to their default values.
@@ -141,6 +144,9 @@ subroutine write_cputime(day, n, CS, nmax, call_end)
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "write_energy: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "write_cputime: Module must be initialized before it is used.")
+
   call SYSTEM_CLOCK(new_cputime, CLOCKS_PER_SEC, MAX_TICKS)
 !   The following lines extract useful information even if the clock has rolled
 ! over, assuming a 32-bit SYSTEM_CLOCK.  With more bits, rollover is essentially
diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 588fa4c75e..94b62f1c20 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -29,6 +29,7 @@ module MOM_MEKE
 
 !> Control structure that contains MEKE parameters and diagnostics handles
 type, public :: MEKE_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   ! Parameters
   real, dimension(:,:), pointer :: equilibrium_value => NULL() !< The equilbrium value
                         !! of MEKE to be calculated at each time step [L2 T-2 ~> m2 s-2]
@@ -180,6 +181,9 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
   if (.not.associated(MEKE)) call MOM_error(FATAL, &
          "MOM_MEKE: MEKE must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_MEKE: Module must be initialized before it is used.")
+
   if ((CS%MEKE_Cd_scale > 0.0) .or. (CS%MEKE_Cb>0.) .or. CS%visc_drag) then
     use_drag_rate = .true.
   else
@@ -1064,6 +1068,8 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   call MOM_mesg("MEKE_init: reading parameters ", 5)
 
   ! Read all relevant parameters and write them to the model log.
diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90
index e0ec7fba63..e3ea65c1a8 100644
--- a/src/parameterizations/lateral/MOM_hor_visc.F90
+++ b/src/parameterizations/lateral/MOM_hor_visc.F90
@@ -30,6 +30,7 @@ module MOM_hor_visc
 
 !> Control structure for horizontal viscosity
 type, public :: hor_visc_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: Laplacian       !< Use a Laplacian horizontal viscosity if true.
   logical :: biharmonic      !< Use a biharmonic horizontal viscosity if true.
   logical :: debug           !< If true, write verbose checksums for debugging purposes.
@@ -421,6 +422,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "MOM_hor_visc: Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_hor_visc: Module must be initialized before it is used.")
+
   if (.not.(CS%Laplacian .or. CS%biharmonic)) return
 
   find_FrictWork = (CS%id_FrictWork > 0)
@@ -1836,6 +1841,9 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp)
     return
   endif
   allocate(CS)
+
+  CS%initialized = .true.
+
   CS%diag => diag
   ! Read parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
index 9306233112..ccd583b165 100644
--- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
+++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
@@ -25,6 +25,7 @@ module MOM_lateral_mixing_coeffs
 
 !> Variable mixing coefficients
 type, public :: VarMix_CS
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: use_variable_mixing  !< If true, use the variable mixing.
   logical :: Resoln_scaling_used  !< If true, a resolution function is used somewhere to scale
                                   !! away one of the viscosities or diffusivities when the
@@ -177,6 +178,10 @@ subroutine calc_depth_function(G, CS)
 
   if (.not. associated(CS)) call MOM_error(FATAL, "calc_depth_function:"// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "calc_depth_function: "// &
+         "Module must be initialized before it is used.")
+
   if (.not. CS%calculate_depth_fns) return
   if (.not. associated(CS%Depth_fn_u)) call MOM_error(FATAL, &
     "calc_depth_function: %Depth_fn_u is not associated with Depth_scaled_KhTh.")
@@ -220,6 +225,10 @@ subroutine calc_resoln_function(h, tv, G, GV, US, CS)
 
   if (.not. associated(CS)) call MOM_error(FATAL, "calc_resoln_function:"// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "calc_resoln_function: "// &
+         "Module must be initialized before it is used.")
+
   if (CS%calculate_cg1) then
     if (.not. associated(CS%cg1)) call MOM_error(FATAL, &
       "calc_resoln_function: %cg1 is not associated with Resoln_scaled_Kh.")
@@ -465,6 +474,9 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC)
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions:"//&
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_lateral_mixing_coeffs.F90, calc_slope_functions: "//&
+         "Module must be initialized before it is used.")
+
   if (CS%calculate_Eady_growth_rate) then
     if (CS%use_simpler_Eady_growth_rate) then
       call find_eta(h, tv, G, GV, US, e, halo_size=2)
@@ -533,6 +545,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C
 
   if (.not. associated(CS)) call MOM_error(FATAL, "calc_slope_function:"// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "calc_Visbeck_coeffs_old: "// &
+         "Module must be initialized before it is used.")
+
   if (.not. CS%calculate_Eady_growth_rate) return
   if (.not. associated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// &
          "%SN_u is not associated with use_variable_mixing.")
@@ -880,6 +896,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop
 
   if (.not. associated(CS)) call MOM_error(FATAL, "calc_slope_function:"// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "calc_slope_functions_using_just_e: "// &
+         "Module must be initialized before it is used.")
+
   if (.not. CS%calculate_Eady_growth_rate) return
   if (.not. associated(CS%SN_u)) call MOM_error(FATAL, "calc_slope_function:"// &
          "%SN_u is not associated with use_variable_mixing.")
@@ -1186,6 +1206,9 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS)
   endif
 
   allocate(CS)
+
+  CS%initialized = .true.
+
   in_use = .false. ! Set to true to avoid deallocating
   CS%diag => diag ! Diagnostics pointer
   CS%calculate_cg1 = .false.
diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
index 0d2062441e..cb1aa16d44 100644
--- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
+++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
@@ -36,6 +36,7 @@ module MOM_mixed_layer_restrat
 
 !> Control structure for mom_mixed_layer_restrat
 type, public :: mixedlayer_restrat_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real    :: ml_restrat_coef       !< A non-dimensional factor by which the instability is enhanced
                                    !! over what would be predicted based on the resolved gradients
                                    !! [nondim].  This increases with grid spacing^2, up to something
@@ -107,6 +108,9 @@ subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, VarMix, G, GV,
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
+         "Module must be initialized before it is used.")
+
   if (GV%nkml>0) then
     call mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
   else
@@ -616,6 +620,10 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_mixedlayer_restrat: "// &
+         "Module must be initialized before it is used.")
+
   if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return
 
   uDml(:)    = 0.0 ; vDml(:) = 0.0
@@ -826,6 +834,8 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
     call MOM_error(FATAL, "mixedlayer_restrat_init called without an associated control structure.")
   endif
 
+  CS%initialized = .true.
+
   ! Nonsense values to cause problems when these parameters are not used
   CS%MLE_MLD_decay_time = -9.e9*US%s_to_T
   CS%MLE_density_diff = -9.e9*US%kg_m3_to_R
diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
index 78425676b1..3f4ca58d0c 100644
--- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90
+++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90
@@ -35,6 +35,7 @@ module MOM_thickness_diffuse
 
 !> Control structure for thickness diffusion
 type, public :: thickness_diffuse_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real    :: Khth                !< Background interface depth diffusivity [L2 T-1 ~> m2 s-1]
   real    :: Khth_Slope_Cff      !< Slope dependence coefficient of Khth [nondim]
   real    :: max_Khth_CFL        !< Maximum value of the diffusive CFL for thickness diffusion
@@ -164,6 +165,9 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_thickness_diffuse: "//&
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_thickness_diffuse: "//&
+         "Module must be initialized before it is used.")
+
   if ((.not.CS%thickness_diffuse) .or. &
        .not.( CS%Khth > 0.0 .or. associated(VarMix) .or. associated(MEKE) ) ) return
 
@@ -1910,6 +1914,8 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
 
   ! Read all relevant parameters and write them to the model log.
diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
index be2dfefe8c..d6c1269614 100644
--- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
+++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90
@@ -30,6 +30,7 @@ module MOM_bulk_mixed_layer
 
 !> The control structure with parameters for the MOM_bulk_mixed_layer module
 type, public :: bulkmixedlayer_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   integer :: nkml            !< The number of layers in the mixed layer.
   integer :: nkbl            !< The number of buffer layers.
   integer :: nsw             !< The number of bands of penetrating shortwave radiation.
@@ -331,6 +332,10 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_mixed_layer: "//&
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_bulk_mixed_layer: "//&
+         "Module must be initialized before it is used.")
+
   if (GV%nkml < 1) return
 
   if (.not. associated(tv%eqn_of_state)) call MOM_error(FATAL, &
@@ -3383,6 +3388,8 @@ subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
   CS%Time => Time
 
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index d3f92e99cc..74c7d48db7 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -88,6 +88,7 @@ module MOM_diabatic_driver
 
 !> Control structure for this module
 type, public :: diabatic_CS ; private
+  logical :: initialized = .false.   !< True if this control structure has been initialized.
 
   logical :: use_legacy_diabatic     !< If true (default), use a legacy version of the diabatic
                                      !! algorithm. This is temporary and is needed to avoid change
@@ -305,6 +306,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
+         "Module must be initialized before it is used.")
+
   if (dt == 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
         "diabatic was called with a zero length timestep.")
   if (dt < 0.0) call MOM_error(FATAL, "MOM_diabatic_driver: "// &
@@ -2908,6 +2913,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
     allocate(CS)
   endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
   CS%Time => Time
 
diff --git a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90 b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
index 2cca587e9e..bd50b4ec0a 100644
--- a/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
+++ b/src/parameterizations/vertical/MOM_diapyc_energy_req.F90
@@ -77,6 +77,9 @@ subroutine diapyc_energy_req_test(h_3d, dt, tv, G, GV, US, CS, Kd_int)
   if (.not. associated(CS)) call MOM_error(FATAL, "diapyc_energy_req_test: "// &
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "diapyc_energy_req_test: "// &
+         "Module must be initialized before it is used.")
+
 !$OMP do
   do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
     if (present(Kd_int) .and. .not.CS%use_test_Kh_profile) then
diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90
index 946a40d39e..8379e1389f 100644
--- a/src/parameterizations/vertical/MOM_energetic_PBL.F90
+++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90
@@ -32,6 +32,7 @@ module MOM_energetic_PBL
 
 !> This control structure holds parameters for the MOM_energetic_PBL module
 type, public :: energetic_PBL_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
 
   !/ Constants
   real    :: VonKar = 0.41   !< The von Karman coefficient.  This should be a runtime parameter,
@@ -348,6 +349,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS
   if (.not. associated(CS)) call MOM_error(FATAL, "energetic_PBL: "//&
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "energetic_PBL: "//&
+         "Module must be initialized before it is used.")
+
   if (.not. associated(tv%eqn_of_state)) call MOM_error(FATAL, &
       "energetic_PBL: Temperature, salinity and an equation of state "//&
       "must now be used.")
@@ -1938,6 +1942,8 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
   CS%Time => Time
 
diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
index 4dc08284af..f6ae73857c 100644
--- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90
+++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90
@@ -27,6 +27,7 @@ module MOM_entrain_diffusive
 
 !> The control structure holding parametes for the MOM_entrain_diffusive module
 type, public :: entrain_diffusive_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: bulkmixedlayer  !< If true, a refined bulk mixed layer is used with
                              !! GV%nk_rho_varies variable density mixed & buffer layers.
   integer :: max_ent_it      !< The maximum number of iterations that may be used to
@@ -210,6 +211,9 @@ subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
   if (.not. associated(CS)) call MOM_error(FATAL, &
          "MOM_entrain_diffusive: Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, &
+         "MOM_entrain_diffusive: Module must be initialized before it is used.")
+
   if (.not.(present(Kd_Lay) .or. present(Kd_int))) call MOM_error(FATAL, &
       "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.")
 
@@ -2087,6 +2091,8 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_re
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   CS%diag => diag
 
   CS%bulkmixedlayer = (GV%nkml > 0)
diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90
index 877f9a0497..802254ee01 100644
--- a/src/parameterizations/vertical/MOM_geothermal.F90
+++ b/src/parameterizations/vertical/MOM_geothermal.F90
@@ -23,6 +23,7 @@ module MOM_geothermal
 
 !> Control structure for geothermal heating
 type, public :: geothermal_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real    :: dRcv_dT_inplace  !< The value of dRcv_dT above which (dRcv_dT is negative) the
                               !! water is heated in place instead of moving upward between
                               !! layers in non-ALE layered mode [R degC-1 ~> kg m-3 degC-1]
@@ -121,6 +122,10 @@ subroutine geothermal_entraining(h, tv, dt, ea, eb, G, GV, US, CS, halo)
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_geothermal: "//&
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//&
+         "Module must be initialized before it is used.")
+
   if (.not.CS%apply_geothermal) return
 
   nkmb      = GV%nk_rho_varies
@@ -397,6 +402,10 @@ subroutine geothermal_in_place(h, tv, dt, G, GV, US, CS, halo)
 
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_geothermal: "//&
          "Module must be initialized before it is used.")
+
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_geothermal: "//&
+         "Module must be initialized before it is used.")
+
   if (.not.CS%apply_geothermal) return
 
   Irho_cp   = 1.0 / (GV%H_to_RZ * tv%C_p)
@@ -518,6 +527,8 @@ subroutine geothermal_init(Time, G, GV, US, param_file, diag, CS, useALEalgorith
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
   CS%Time => Time
 
diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90
index df24d3f4e9..300cdcbe1e 100644
--- a/src/parameterizations/vertical/MOM_internal_tide_input.F90
+++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90
@@ -34,6 +34,7 @@ module MOM_int_tide_input
 
 !> This control structure holds parameters that regulate internal tide energy inputs.
 type, public :: int_tide_input_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: debug      !< If true, write verbose checksums for debugging.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                         !! regulate the timing of diagnostic output.
@@ -111,6 +112,9 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
   if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//&
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"set_diffusivity: "//&
+         "Module must be initialized before it is used.")
+
   use_EOS = associated(tv%eqn_of_state)
 
   ! Smooth the properties through massless layers.
@@ -330,6 +334,8 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
   allocate(CS)
   allocate(itide)
 
+  CS%initialized = .true.
+
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
diff --git a/src/parameterizations/vertical/MOM_regularize_layers.F90 b/src/parameterizations/vertical/MOM_regularize_layers.F90
index af92e522a2..79d6d064c7 100644
--- a/src/parameterizations/vertical/MOM_regularize_layers.F90
+++ b/src/parameterizations/vertical/MOM_regularize_layers.F90
@@ -23,6 +23,7 @@ module MOM_regularize_layers
 
 !> This control structure holds parameters used by the MOM_regularize_layers module
 type, public :: regularize_layers_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: regularize_surface_layers !< If true, vertically restructure the
                              !! near-surface layers when they have too much
                              !! lateral variations to allow for sensible lateral
@@ -96,6 +97,9 @@ subroutine regularize_layers(h, tv, dt, ea, eb, G, GV, US, CS)
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_regularize_layers: "//&
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_regularize_layers: "//&
+         "Module must be initialized before it is used.")
+
   if (CS%regularize_surface_layers) then
     call pass_var(h, G%Domain, clock=id_clock_pass)
     call regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS)
@@ -197,6 +201,9 @@ subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS)
   if (.not. associated(CS)) call MOM_error(FATAL, "MOM_regularize_layers: "//&
          "Module must be initialized before it is used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_regularize_layers: "//&
+         "Module must be initialized before it is used.")
+
   if (GV%nkml<1) return
   nkmb = GV%nk_rho_varies ; nkml = GV%nkml
   if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, &
@@ -735,6 +742,8 @@ subroutine regularize_layers_init(Time, G, GV, param_file, diag, CS)
     return
   else ; allocate(CS) ; endif
 
+  CS%initialized = .true.
+
   CS%diag => diag
   CS%Time => Time
 
diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90
index b1a4d1433d..89a426d0cd 100644
--- a/src/parameterizations/vertical/MOM_set_diffusivity.F90
+++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90
@@ -56,6 +56,7 @@ module MOM_set_diffusivity
 
 !> This control structure contains parameters for MOM_set_diffusivity.
 type, public :: set_diffusivity_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: debug           !< If true, write verbose checksums for debugging.
 
   logical :: bulkmixedlayer  !< If true, a refined bulk mixed layer is used with
@@ -282,6 +283,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
   if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//&
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"set_diffusivity: "//&
+         "Module must be initialized before it is used.")
+
   if (CS%answers_2018) then
     ! These hard-coded dimensional parameters are being replaced.
     kappa_dt_fill = US%m_to_Z**2 * 1.e-3 * 7200.
@@ -1707,6 +1711,9 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC)
   if (.not.associated(CS)) call MOM_error(FATAL,"set_BBL_TKE: "//&
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"set_BBL_TKE: "//&
+         "Module must be initialized before it is used.")
+
   if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0)) then
     if (associated(visc%ustar_BBL)) then
       do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo
@@ -1986,6 +1993,8 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   is  = G%isc ; ie  = G%iec ; js  = G%jsc ; je  = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
 
diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90
index 1cf3b5ddc9..ca7964c43c 100644
--- a/src/parameterizations/vertical/MOM_set_viscosity.F90
+++ b/src/parameterizations/vertical/MOM_set_viscosity.F90
@@ -43,6 +43,7 @@ module MOM_set_visc
 
 !> Control structure for MOM_set_visc
 type, public :: set_visc_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real    :: Hbbl           !< The static bottom boundary layer thickness [H ~> m or kg m-2].
                             !! Runtime parameter `HBBL`.
   real    :: cdrag          !< The quadratic drag coefficient.
@@ -286,6 +287,10 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//&
          "Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//&
+         "Module must be initialized before it is used.")
+
   if (.not.CS%bottomdraglaw) return
 
   if (CS%debug) then
@@ -1249,6 +1254,10 @@ subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
 
   if (.not.associated(CS)) call MOM_error(FATAL,"MOM_set_viscosity(visc_ML): "//&
          "Module must be initialized before it is used.")
+
+  if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(visc_ML): "//&
+         "Module must be initialized before it is used.")
+
   if (.not.(CS%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
             associated(forces%frac_shelf_v)) ) return
 
@@ -1939,6 +1948,8 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   CS%OBC => OBC
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90
index d0d3943a26..f36af41149 100644
--- a/src/parameterizations/vertical/MOM_vert_friction.F90
+++ b/src/parameterizations/vertical/MOM_vert_friction.F90
@@ -37,6 +37,7 @@ module MOM_vert_friction
 
 !> The control structure with parameters and memory for the MOM_vert_friction module
 type, public :: vertvisc_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real    :: Hmix            !< The mixed layer thickness in thickness units [H ~> m or kg m-2].
   real    :: Hmix_stress     !< The mixed layer thickness over which the wind
                              !! stress is applied with direct_stress [H ~> m or kg m-2].
@@ -236,6 +237,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
   if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
+         "Module must be initialized before it is used.")
+
   if (CS%direct_stress) then
     Hmix = CS%Hmix_stress
     I_Hmix = 1.0 / Hmix
@@ -648,6 +652,9 @@ subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
   if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(remant): "// &
+         "Module must be initialized before it is used.")
+
   dt_Z_to_H = dt*GV%Z_to_H
 
   do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo
@@ -800,6 +807,9 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, US, CS, OBC)
   if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(coef): "// &
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(coef): "// &
+         "Module must be initialized before it is used.")
+
   h_neglect = GV%H_subroundoff
   a_cpl_max = 1.0e37 * US%m_to_Z * US%T_to_s
   I_Hbbl(:) = 1.0 / (CS%Hbbl + h_neglect)
@@ -1704,6 +1714,8 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   if (GV%Boussinesq) then; thickness_units = "m"
   else; thickness_units = "kg m-2"; endif
 
diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90
index 92570e3caa..0308a3b008 100644
--- a/src/user/user_change_diffusivity.F90
+++ b/src/user/user_change_diffusivity.F90
@@ -25,6 +25,7 @@ module user_change_diffusivity
 
 !> Control structure for user_change_diffusivity
 type, public :: user_change_diff_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   real :: Kd_add        !< The scale of a diffusivity that is added everywhere
                         !! without any filtering or scaling [Z2 T-1 ~> m2 s-1].
   real :: lat_range(4)  !< 4 values that define the latitude range over which
@@ -86,6 +87,9 @@ subroutine user_change_diff(h, tv, G, GV, US, CS, Kd_lay, Kd_int, T_f, S_f, Kd_i
   if (.not.associated(CS)) call MOM_error(FATAL,"user_set_diffusivity: "//&
          "Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL,"user_set_diffusivity: "//&
+         "Module must be initialized before it is used.")
+
   use_EOS = associated(tv%eqn_of_state)
   if (.not.use_EOS) return
   store_Kd_add = .false.
@@ -214,6 +218,8 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS)
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
 
   CS%diag => diag