diff --git a/opm/models/blackoil/blackoilenergymodules.hh b/opm/models/blackoil/blackoilenergymodules.hh
index e8007d3d4a..c09fe92439 100644
--- a/opm/models/blackoil/blackoilenergymodules.hh
+++ b/opm/models/blackoil/blackoilenergymodules.hh
@@ -54,13 +54,16 @@
 
 namespace Opm {
 
+template <class TypeTag, EnergyModules activeModule>
+class BlackOilEnergyModule;
+
 /*!
  * \ingroup BlackOil
  * \brief Contains the high level supplements required to extend the black oil
  *        model by energy.
  */
-template <class TypeTag, EnergyModules activeModule = getPropValue<TypeTag, Properties::EnergyModuleType>()>
-class BlackOilEnergyModule
+template <class TypeTag>
+class BlackOilEnergyModule<TypeTag, EnergyModules::FullyImplicitThermal>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@@ -77,7 +80,7 @@ class BlackOilEnergyModule
     static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
     static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
 
-    static constexpr unsigned enableFullyImplicitThermal = (activeModule == EnergyModules::FullyImplicitThermal);
+    static constexpr unsigned enableFullyImplicitThermal = true;
     static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
     static constexpr unsigned numPhases = FluidSystem::numPhases;
 
@@ -89,9 +92,7 @@ public:
      */
     static void registerParameters()
     {
-        if constexpr (enableFullyImplicitThermal) {
-            VtkBlackOilEnergyModule<TypeTag>::registerParameters();
-        }
+        VtkBlackOilEnergyModule<TypeTag>::registerParameters();
     }
 
     /*!
@@ -100,19 +101,12 @@ public:
     static void registerOutputModules(Model& model,
                                       Simulator& simulator)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            model.addOutputModule(std::make_unique<VtkBlackOilEnergyModule<TypeTag>>(simulator));
-        }
+        model.addOutputModule(std::make_unique<VtkBlackOilEnergyModule<TypeTag>>(simulator));
     }
 
     static OPM_HOST_DEVICE bool primaryVarApplies(unsigned pvIdx)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            return pvIdx == temperatureIdx;
-        }
-        else {
-            return false;
-        }
+        return pvIdx == temperatureIdx;
     }
 
     static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
@@ -132,12 +126,7 @@ public:
 
     static OPM_HOST_DEVICE bool eqApplies(unsigned eqIdx)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            return eqIdx == contiEnergyEqIdx;
-        }
-        else {
-            return false;
-        }
+        return eqIdx == contiEnergyEqIdx;
     }
 
     static std::string eqName([[maybe_unused]] unsigned eqIdx)
@@ -160,73 +149,66 @@ public:
                                            const IntensiveQuantities& intQuants)
     {
         using LhsEval = typename StorageType::value_type;
+        const FluidSystem& fsys = intQuants.getFluidSystem();
 
-        if constexpr (enableFullyImplicitThermal) {
-            const FluidSystem& fsys = intQuants.getFluidSystem();
+        const auto& poro = decay<LhsEval>(intQuants.porosity());
 
-            const auto& poro = decay<LhsEval>(intQuants.porosity());
-
-            // accumulate the internal energy of the fluids
-            const auto& fs = intQuants.fluidState();
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-                if (!fsys.phaseIsActive(phaseIdx)) {
-                    continue;
-                }
-
-                const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
-                const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
-                const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
-
-                storage[contiEnergyEqIdx] += poro*S*u*rho;
+        // accumulate the internal energy of the fluids
+        const auto& fs = intQuants.fluidState();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!fsys.phaseIsActive(phaseIdx)) {
+                continue;
             }
 
-            // add the internal energy of the rock
-            const Scalar rockFraction = intQuants.rockFraction();
-            const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
-            storage[contiEnergyEqIdx] += rockFraction * uRock;
-            storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+            const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
+            const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
+            const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
+
+            storage[contiEnergyEqIdx] += poro*S*u*rho;
         }
+
+        // add the internal energy of the rock
+        const Scalar rockFraction = intQuants.rockFraction();
+        const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
+        storage[contiEnergyEqIdx] += rockFraction * uRock;
+        storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
-    OPM_HOST_DEVICE static void computeFlux([[maybe_unused]] RateVector& flux,
-                                            [[maybe_unused]] const ElementContext& elemCtx,
-                                            [[maybe_unused]] unsigned scvfIdx,
-                                            [[maybe_unused]] unsigned timeIdx)
+    OPM_HOST_DEVICE static void computeFlux(RateVector& flux,
+                                            const ElementContext& elemCtx,
+                                            unsigned scvfIdx,
+                                            unsigned timeIdx)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            flux[contiEnergyEqIdx] = 0.0;
+        flux[contiEnergyEqIdx] = 0.0;
 
-            const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
-            const unsigned focusIdx = elemCtx.focusDofIndex();
-            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                if (!FluidSystem::phaseIsActive(phaseIdx)) {
-                    continue;
-                }
-
-                const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
-                if (upIdx == focusIdx) {
-                    addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
-                }
-                else {
-                    addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
-                }
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        const unsigned focusIdx = elemCtx.focusDofIndex();
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
             }
 
-            // diffusive energy flux
-            flux[contiEnergyEqIdx] += extQuants.energyFlux();
-            flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+            const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
+            if (upIdx == focusIdx) {
+                addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
+            }
+            else {
+                addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
+            }
         }
+
+        // diffusive energy flux
+        flux[contiEnergyEqIdx] += extQuants.energyFlux();
+        flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     template<class RateVectorT>
     OPM_HOST_DEVICE static void addHeatFlux(RateVectorT& flux,
                                             const Evaluation& heatFlux)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            // diffusive energy flux
-            flux[contiEnergyEqIdx] += heatFlux;
-            flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
-        }
+        // diffusive energy flux
+        flux[contiEnergyEqIdx] += heatFlux;
+        flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     template <class UpEval, class RateVectorT, class Eval, class FluidState>
@@ -262,9 +244,7 @@ public:
     OPM_HOST_DEVICE static void addToEnthalpyRate(RateVector& flux,
                                                   const Evaluation& hRate)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            flux[contiEnergyEqIdx] += hRate;
-        }
+        flux[contiEnergyEqIdx] += hRate;
     }
 
     /*!
@@ -274,9 +254,7 @@ public:
     OPM_HOST_DEVICE static void assignPrimaryVars(PrimaryVariables& priVars,
                                                   const FluidState& fluidState)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            priVars[temperatureIdx] = getValue(fluidState.temperature(/*phaseIdx=*/0));
-        }
+        priVars[temperatureIdx] = getValue(fluidState.temperature(/*phaseIdx=*/0));
     }
 
     /*!
@@ -286,10 +264,8 @@ public:
                                                   const PrimaryVariables& oldPv,
                                                   const EqVector& delta)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            // do a plain unchopped Newton update
-            newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
-        }
+        // do a plain unchopped Newton update
+        newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
     }
 
     /*!
@@ -316,26 +292,22 @@ public:
     template <class DofEntity>
     static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            const unsigned dofIdx = model.dofMapper().index(dof);
-            const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
-            outstream << priVars[temperatureIdx];
-        }
+        const unsigned dofIdx = model.dofMapper().index(dof);
+        const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
+        outstream << priVars[temperatureIdx];
     }
 
     template <class DofEntity>
     static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
     {
-        if constexpr (enableFullyImplicitThermal) {
-            const unsigned dofIdx = model.dofMapper().index(dof);
-            PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
-            PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
+        const unsigned dofIdx = model.dofMapper().index(dof);
+        PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
+        PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
 
-            instream >> priVars0[temperatureIdx];
+        instream >> priVars0[temperatureIdx];
 
-            // set the primary variables for the beginning of the current time step.
-            priVars1 = priVars0[temperatureIdx];
-        }
+        // set the primary variables for the beginning of the current time step.
+        priVars1 = priVars0[temperatureIdx];
     }
 };
 
diff --git a/opm/models/blackoil/blackoillocalresidual.hh b/opm/models/blackoil/blackoillocalresidual.hh
index 17f7feb876..5e15f025ed 100644
--- a/opm/models/blackoil/blackoillocalresidual.hh
+++ b/opm/models/blackoil/blackoillocalresidual.hh
@@ -83,8 +83,8 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
     static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
     static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
     static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
-    static constexpr bool enableFullyImplicitThermal = 
-        getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal;
+    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
+    static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
     static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
     static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
 
@@ -94,7 +94,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
     using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
     using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
-    using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
     using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
     using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
     using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
@@ -191,7 +191,9 @@ public:
         }
 
         // deal with energy (if present)
-        EnergyModule::addStorage(storage, intQuants);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::addStorage(storage, intQuants);
+        }
 
         // deal with foam (if present)
         if constexpr (enableFoam) {
@@ -255,7 +257,9 @@ public:
         }
 
         // deal with energy (if present)
-        EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+        }
 
         // deal with foam (if present)
         if constexpr (enableFoam) {
diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh
index 612954ecdb..5d3de0b6f0 100644
--- a/opm/models/blackoil/blackoillocalresidualtpfa.hh
+++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh
@@ -101,9 +101,10 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
     static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
     static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
-    static constexpr bool enableFullyImplicitThermal
-        = (getPropValue<TypeTag, Properties::EnergyModuleType>()
-           == EnergyModules::FullyImplicitThermal);
+    static constexpr EnergyModules energyModuleType =
+        getPropValue<TypeTag, Properties::EnergyModuleType>();
+    static constexpr bool enableFullyImplicitThermal =
+        energyModuleType == EnergyModules::FullyImplicitThermal;
     static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
     static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
     static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
@@ -121,7 +122,7 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
     using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
-    using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
     using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
     using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
     using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
@@ -233,7 +234,9 @@ public:
         }
 
         // deal with energy (if present)
-        EnergyModule::addStorage(storage, intQuants);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::addStorage(storage, intQuants);
+        }
 
         // deal with foam (if present)
         if constexpr (enableFoam) {
diff --git a/opm/models/blackoil/blackoilmodel.hh b/opm/models/blackoil/blackoilmodel.hh
index ec5111da00..486064501d 100644
--- a/opm/models/blackoil/blackoilmodel.hh
+++ b/opm/models/blackoil/blackoilmodel.hh
@@ -358,12 +358,14 @@ private:
     static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
     static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
     static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
+    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
+    static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
     static constexpr bool waterEnabled = Indices::waterEnabled;
 
     using BioeffectsModule = BlackOilBioeffectsModule<TypeTag, enableBioeffects>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
     using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
-    using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
     using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
     using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
     using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
@@ -393,7 +395,9 @@ public:
         if constexpr (enablePolymer) {
             PolymerModule::registerParameters();
         }
-        EnergyModule::registerParameters();
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::registerParameters();
+        }
         if constexpr (enableDiffusion) {
             DiffusionModule::registerParameters();
         }
@@ -448,12 +452,13 @@ public:
             }
         }
 
-        if (EnergyModule::primaryVarApplies(pvIdx)) {
-            return EnergyModule::primaryVarName(pvIdx);
-        }
-        else {
-            throw std::logic_error("Invalid primary variable index");
+        if constexpr (enableFullyImplicitThermal) {
+            if (EnergyModule::primaryVarApplies(pvIdx)) {
+                return EnergyModule::primaryVarName(pvIdx);
+            }
         }
+
+        throw std::logic_error("Invalid primary variable index");
     }
 
     /*!
@@ -485,12 +490,13 @@ public:
             }
         }
 
-        if (EnergyModule::eqApplies(eqIdx)) {
-            return EnergyModule::eqName(eqIdx);
-        }
-        else {
-            throw std::logic_error("Invalid equation index");
+        if constexpr (enableFullyImplicitThermal) {
+            if (EnergyModule::eqApplies(eqIdx)) {
+                return EnergyModule::eqName(eqIdx);
+            }
         }
+
+        throw std::logic_error("Invalid equation index");
     }
 
     /*!
@@ -537,8 +543,10 @@ public:
         }
 
         // deal with primary variables stemming from the energy module
-        if (EnergyModule::primaryVarApplies(pvIdx)) {
-            return EnergyModule::primaryVarWeight(pvIdx);
+        if constexpr (enableFullyImplicitThermal) {
+            if (EnergyModule::primaryVarApplies(pvIdx)) {
+                return EnergyModule::primaryVarWeight(pvIdx);
+            }
         }
 
         // if the primary variable is either the gas saturation, Rs or Rv
@@ -612,7 +620,9 @@ public:
         if constexpr (enablePolymer) {
             PolymerModule::serializeEntity(asImp_(), outstream, dof);
         }
-        EnergyModule::serializeEntity(asImp_(), outstream, dof);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::serializeEntity(asImp_(), outstream, dof);
+        }
     }
 
     /*!
@@ -664,7 +674,9 @@ public:
         if constexpr (enablePolymer) {
             PolymerModule::deserializeEntity(asImp_(), instream, dof);
         }
-        EnergyModule::deserializeEntity(asImp_(), instream, dof);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::deserializeEntity(asImp_(), instream, dof);
+        }
 
         using PVM_G = typename PrimaryVariables::GasMeaning;
         using PVM_W = typename PrimaryVariables::WaterMeaning;
@@ -731,7 +743,9 @@ protected:
         if constexpr (enablePolymer) {
             PolymerModule::registerOutputModules(asImp_(), this->simulator_);
         }
-        EnergyModule::registerOutputModules(asImp_(), this->simulator_);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::registerOutputModules(asImp_(), this->simulator_);
+        }
         if constexpr (enableBioeffects) {
             BioeffectsModule::registerOutputModules(asImp_(), this->simulator_);
         }
diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh
index 5525dccaff..40df6c452d 100644
--- a/opm/models/blackoil/blackoilprimaryvariables.hh
+++ b/opm/models/blackoil/blackoilprimaryvariables.hh
@@ -102,16 +102,18 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag, VectorTy
     // component indices from the fluid system
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
 
+    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
+
     static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
     static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
     static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
     static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
+    static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
     static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
     static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
 
     enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
     enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
-    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     enum { enableMICP = Indices::enableMICP };
     enum { gasCompIdx = FluidSystem::gasCompIdx };
     enum { waterCompIdx = FluidSystem::waterCompIdx };
@@ -339,7 +341,9 @@ public:
         const bool precipitatedSaltPresent = enableSaltPrecipitation ? saltSaturation > 0.0 : false;
         const bool oneActivePhases = fluidState.fluidSystem().numActivePhases() == 1;
         // deal with the primary variables for the energy extension
-        EnergyModule::assignPrimaryVars(*this, fluidState);
+        if constexpr (enableFullyImplicitThermal) {
+            EnergyModule::assignPrimaryVars(*this, fluidState);
+        }
 
         // Determine the meaning of the pressure primary variables
         // Depending on the phases present, this variable is either interpreted as the
diff --git a/opm/simulators/flow/TemperatureModel.hpp b/opm/simulators/flow/TemperatureModel.hpp
index 0155ebcf91..8a5d3ad344 100644
--- a/opm/simulators/flow/TemperatureModel.hpp
+++ b/opm/simulators/flow/TemperatureModel.hpp
@@ -196,7 +196,8 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
     using Indices = GetPropType<TypeTag, Properties::Indices>;
     using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
-    using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
+    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
     using IndexTraits = typename FluidSystem::IndexTraitsType;
     using WellStateType = WellState<Scalar, IndexTraits>;
 
@@ -554,7 +555,7 @@ protected:
         flux *= scalingFactor_;
     }
 
-    template < class ResidualNBInfo>
+    template <class ResidualNBInfo>
     void computeHeatFluxTerm(const IntensiveQuantitiesTemp& intQuantsIn,
                              const IntensiveQuantitiesTemp& intQuantsEx,
                              const ResidualNBInfo& res_nbinfo,
@@ -562,17 +563,18 @@ protected:
     {
         short interiorDofIdx = 0; // NB
         short exteriorDofIdx = 1; // NB
-        EnergyModule::ExtensiveQuantities::updateEnergy(heatFlux,
-                                                        interiorDofIdx, // focusDofIndex,
-                                                        interiorDofIdx,
-                                                        exteriorDofIdx,
-                                                        intQuantsIn,
-                                                        intQuantsEx,
-                                                        intQuantsIn.fluidStateTemp(),
-                                                        intQuantsEx.fluidStateTemp(),
-                                                        res_nbinfo.inAlpha,
-                                                        res_nbinfo.outAlpha,
-                                                        res_nbinfo.faceArea);
+        BlackOilEnergyExtensiveQuantities<TypeTag, energyModuleType>::
+            updateEnergy(heatFlux,
+                         interiorDofIdx, // focusDofIndex,
+                         interiorDofIdx,
+                         exteriorDofIdx,
+                         intQuantsIn,
+                         intQuantsEx,
+                         intQuantsIn.fluidStateTemp(),
+                         intQuantsEx.fluidStateTemp(),
+                         res_nbinfo.inAlpha,
+                         res_nbinfo.outAlpha,
+                         res_nbinfo.faceArea);
         heatFlux *= scalingFactor_*res_nbinfo.faceArea;
     }