diff --git a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp
index 61aff7740e..b8e303268d 100644
--- a/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp
+++ b/src/global_legacy_module/4C_global_legacy_module_validmaterials.cpp
@@ -2873,23 +2873,45 @@ std::unordered_map<Core::Materials::MaterialType, Core::IO::InputSpec> Global::v
 
   /*----------------------------------------------------------------------*/
   {
+    using namespace Core::IO::InputSpecBuilders::Validators;
+
     known_materials[Core::Materials::mvl_reformulated_Johnson_Cook] =
         group("MAT_ViscoplasticLawReformulatedJohnsonCook",
             {
                 parameter<double>("STRAIN_RATE_PREFAC",
-                    {.description = "reference plastic strain rate $\\dot{P}_0$ "}),
+                    {.description = "reference plastic strain rate $\\dot{P}_0$",
+                        .validator = positive<double>()}),
                 parameter<double>("STRAIN_RATE_EXP_FAC",
-                    {.description = "exponential factor of plastic strain rate $C$"}),
+                    {.description = "exponential factor of plastic strain rate $C$",
+                        .validator = positive<double>()}),
                 parameter<double>("INIT_YIELD_STRENGTH",
-                    {.description = "initial yield strength of the material $A_0$"}),
+                    {.description = "initial yield strength of the material $A_0$",
+                        .validator = positive<double>()}),
                 parameter<double>("ISOTROP_HARDEN_PREFAC",
-                    {.description = "prefactor of the isotropic hardening stress $B_0$"}),
+                    {.description =
+                            "prefactor of the isotropic hardening stress / hardening modulus $B_0$",
+                        .validator = positive_or_zero<double>()}),
                 parameter<double>("ISOTROP_HARDEN_EXP",
-                    {.description = "exponent of the isotropic hardening stress $n$"}),
+                    {.description = "exponent of the isotropic hardening stress $n$",
+                        .validator = positive_or_zero<double>()}),
+                parameter<double>("REF_TEMPERATURE",
+                    {.description = "reference temperature $T_0$ for evaluating the "
+                                    "yield strength $A_0$ and the hardening "
+                                    "modulus $B_0$",
+                        .validator = positive<double>()}),
+                parameter<double>(
+                    "MELT_TEMPERATURE", {.description = "melting temperature $T_{\\mathrm{melt}}$",
+                                            .validator = positive<double>()}),
+                parameter<double>("TEMPERATURE_SENS", {.description = "temperature sensitivity $m$",
+                                                          .validator = positive_or_zero<double>()}),
+
             },
             {.description = "Reformulation of the Johnson-Cook viscoplastic law (comprising flow "
-                            "rule $\\dot{P} = \\dot{P}_0 \\exp \\left( \\frac{ \\Sigma_{eq}}{C "
-                            "\\Sigma_y} - \\frac{1}{C} \\right) - \\dot{P}_0$ and hardening law), "
+                            "rule $\\dot{P} = \\dot{P}_0 \\exp \\left( \\frac{ \\sigma_{eq}}{C "
+                            "\\sigma_{\\mathrm{Y}}} - \\frac{1}{C} \\right) - \\dot{P}_0$ and "
+                            "hardening law $\\sigma_{\\mathrm{Y}} = (A_0 + "
+                            "B_0 \\cdot P^{n}) \\cdot (1 - \\frac{T^m - "
+                            "T_{0}^m}{T_{\\mathrm{melt}}^m - T_{0}^m})$), "
                             "as shown in Mareau et al. (Mechanics of Materials 143, 2020)"});
   }
 
diff --git a/src/mat/4C_mat_inelastic_defgrad_factors.cpp b/src/mat/4C_mat_inelastic_defgrad_factors.cpp
index 5c212f49e7..451d782754 100644
--- a/src/mat/4C_mat_inelastic_defgrad_factors.cpp
+++ b/src/mat/4C_mat_inelastic_defgrad_factors.cpp
@@ -1643,6 +1643,9 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::pre_evaluate(
     const Teuchos::ParameterList& params, const EvaluationContext& context, const int gp,
     const int eleGID)
 {
+  // save parameter list
+  params_ = params;
+
   // set Gauss Point
   gp_ = gp;
 
@@ -1662,7 +1665,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::pre_evaluate(
   // call pre_evaluate method of the time step quantities
   time_step_quantities_.pre_evaluate(gp);
   // call preevaluate method of the viscoplastic law
-  viscoplastic_law_->pre_evaluate(gp);
+  viscoplastic_law_->pre_evaluate(params, gp);
 }
 
 /*--------------------------------------------------------------------*
@@ -2231,7 +2234,7 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_state_quantity_deriv
   dNpdMe_sym_dev = Core::LinAlg::Voigt::modify_voigt_representation(temp6x6, 1.0, 2.0);
 
   // compute the relevant derivatives of the plastic strain rate
-  Core::LinAlg::Matrix<2, 1> evoEqFunctionDers =
+  InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs evoEqFunctionDers =
       viscoplastic_law_->evaluate_derivatives_of_plastic_strain_rate(equiv_stress, plastic_strain,
           dt, parameter()->max_plastic_strain_deriv_incr(), err_status);
 
@@ -2249,8 +2252,8 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_state_quantity_deriv
   }
 
   // compute derivatives of the plastic strain rate
-  state_quantity_derivatives.curr_dpsr_dequiv_stress = evoEqFunctionDers(0);
-  state_quantity_derivatives.curr_dpsr_depsp = evoEqFunctionDers(1);
+  state_quantity_derivatives.curr_dpsr_dequiv_stress = evoEqFunctionDers.deriv_equiv_stress;
+  state_quantity_derivatives.curr_dpsr_depsp = evoEqFunctionDers.deriv_plastic_strain;
 
 
   if (eval_type == ViscoplastUtils::StateQuantityDerivEvalType::PlasticStrainRateDerivsOnly)
@@ -2427,7 +2430,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_additional_cmat
             time_step_quantities_.current_plastic_strain[gp_]);
 
     Core::LinAlg::Matrix<10, 10> jacMat(Core::LinAlg::Initialization::zero);
-    viscoplastic_law_->pre_evaluate(gp_);  // set last_substep <- last_
+    viscoplastic_law_->pre_evaluate(params_, gp_);  // set last_substep <- last_
     jacMat = calculate_jacobian(CredM, current_sol,
         time_step_quantities_.last_plastic_defgrad_inverse[gp_],
         time_step_quantities_.last_plastic_strain[gp_], time_step_tracker_.dt, err_status);
@@ -2565,6 +2568,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_inverse_inelast
     {
       time_step_quantities_.current_plastic_defgrad_inverse[gp_] = iFinM;
       time_step_quantities_.current_plastic_strain[gp_] = plastic_strain_pred;
+      time_step_quantities_.current_equiv_stress[gp_] = state_quantities_.curr_equiv_stress;
       time_step_quantities_.current_rightCG[gp_] = CredM;
       time_step_quantities_.current_defgrad[gp_] = FredM;
     }
@@ -2596,6 +2600,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_inverse_inelast
     {
       time_step_quantities_.current_plastic_defgrad_inverse[gp_] = iFinM;
       time_step_quantities_.current_plastic_strain[gp_] = sol(9);
+      time_step_quantities_.current_equiv_stress[gp_] = state_quantities_.curr_equiv_stress;
       time_step_quantities_.current_rightCG[gp_] = CredM;
       time_step_quantities_.current_defgrad[gp_] = FredM;
     }
@@ -2714,7 +2719,7 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::calculate_local_newton_loop_r
   Core::LinAlg::Matrix<3, 3> resFM(Core::LinAlg::Initialization::zero);
   double resepsp = 0.0;
 
-  // compute residuals (standard substepping)
+  // compute residuals (standard time integration)
   if (parameter()->timint_type() == ViscoplastUtils::TimIntType::standard)
   {
     // calculate residual of the equation for inelastic defgrad
@@ -2733,11 +2738,31 @@ Mat::InelasticDefgradTransvIsotropElastViscoplast::calculate_local_newton_loop_r
     last_FinM.invert(last_iFinM);
     Core::LinAlg::Matrix<3, 3> T(Core::LinAlg::Initialization::zero);
     T.multiply_nn(1.0, last_FinM, iFinM, 0.0);
-    Core::LinAlg::MatrixFunctErrorType log_err_status{
-        Core::LinAlg::MatrixFunctErrorType::no_errors};
-    Core::LinAlg::Matrix<3, 3> logT = Core::LinAlg::matrix_log(T, log_err_status);
-    FOUR_C_ASSERT_ALWAYS(log_err_status == Core::LinAlg::MatrixFunctErrorType::no_errors,
-        "Matrix logarithm evaluation failed!");
+    Core::LinAlg::MatrixFunctErrorType log_err_status =
+        Core::LinAlg::MatrixFunctErrorType::no_errors;
+    Core::LinAlg::Matrix<3, 3> logT{Core::LinAlg::Initialization::zero};
+    if (parameter()->mat_log_calc_method() ==
+        Core::LinAlg::MatrixLogCalcMethod::inv_scal_square)  // evaluation using the inverse
+                                                             // scaling and squaring
+                                                             // method
+    {
+      // when computing the matrix logarithm with the inverse scaling
+      // and squaring, we also save the resulting Pade
+      // order via the dedicated pointer. This will be helpful when we
+      // compute the derivative - we want consistent Pade orders for the
+      // evaluations of functions and their derivatives.
+      logT = Core::LinAlg::matrix_log(
+          T, log_err_status, matrix_exp_log_utils_.pade_order, parameter()->mat_log_calc_method());
+    }
+    else  // evaluation using other provided methods
+    {
+      logT = Core::LinAlg::matrix_log(T, log_err_status, parameter()->mat_log_calc_method());
+    }
+    if (log_err_status != Core::LinAlg::MatrixFunctErrorType::no_errors)
+    {
+      err_status = ViscoplastUtils::ErrorType::failed_matrix_log_evaluation;
+      return Core::LinAlg::Matrix<10, 1>{Core::LinAlg::Initialization::zero};
+    }
 
     // calculate residual of the equation for inelastic defgrad
     resFM.update(1.0, logT, dt, state_quantities_.curr_lpM, 0.0);
@@ -3313,5 +3338,88 @@ std::string Mat::InelasticDefgradTransvIsotropElastViscoplast::get_error_info(
 }
 
 
+/*--------------------------------------------------------------------*
+ *--------------------------------------------------------------------*/
+void Mat::InelasticDefgradTransvIsotropElastViscoplast::register_output_data_names(
+    std::unordered_map<std::string, int>& names_and_size) const
+{
+  names_and_size["inverse_plastic_defgrad"] = 9;
+  names_and_size["plastic_strain"] = 1;
+  names_and_size["equiv_stress"] = 1;
+  names_and_size["defgrad"] = 9;
+  names_and_size["rightCG"] = 9;
+  viscoplastic_law_->register_output_data_names(names_and_size);
+}
+
+/*--------------------------------------------------------------------*
+ *--------------------------------------------------------------------*/
+bool Mat::InelasticDefgradTransvIsotropElastViscoplast::evaluate_output_data(
+    const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+{
+  // auxiliaries
+  Core::LinAlg::Matrix<9, 1> temp9x1{Core::LinAlg::Initialization::zero};
+
+  if (name == "inverse_plastic_defgrad")
+  {
+    for (int gp = 0;
+        gp < static_cast<int>(time_step_quantities_.current_plastic_defgrad_inverse.size()); ++gp)
+    {
+      Core::LinAlg::Voigt::matrix_3x3_to_9x1(
+          time_step_quantities_.current_plastic_defgrad_inverse[gp], temp9x1);
+
+      for (int col = 0; col < 9; ++col)
+      {
+        data(gp, col) = temp9x1(col);
+      }
+    }
+    return true;
+  }
+  else if (name == "plastic_strain")
+  {
+    for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_plastic_strain.size());
+        ++gp)
+    {
+      data(gp, 0) = time_step_quantities_.current_plastic_strain[gp];
+    }
+    return true;
+  }
+  else if (name == "equiv_stress")
+  {
+    for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_equiv_stress.size()); ++gp)
+    {
+      data(gp, 0) = time_step_quantities_.current_equiv_stress[gp];
+    }
+    return true;
+  }
+  else if (name == "defgrad")
+  {
+    for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_defgrad.size()); ++gp)
+    {
+      Core::LinAlg::Voigt::matrix_3x3_to_9x1(time_step_quantities_.current_defgrad[gp], temp9x1);
+
+      for (int col = 0; col < 9; ++col)
+      {
+        data(gp, col) = temp9x1(col);
+      }
+    }
+    return true;
+  }
+  else if (name == "rightCG")
+  {
+    for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_rightCG.size()); ++gp)
+    {
+      Core::LinAlg::Voigt::matrix_3x3_to_9x1(time_step_quantities_.current_rightCG[gp], temp9x1);
+
+
+      for (int col = 0; col < 9; ++col)
+      {
+        data(gp, col) = temp9x1(col);
+      }
+    }
+    return true;
+  }
+
+  return viscoplastic_law_->evaluate_output_data(name, data);
+}
 
 FOUR_C_NAMESPACE_CLOSE
diff --git a/src/mat/4C_mat_inelastic_defgrad_factors.hpp b/src/mat/4C_mat_inelastic_defgrad_factors.hpp
index 3eb3b3f684..f656bf96b3 100644
--- a/src/mat/4C_mat_inelastic_defgrad_factors.hpp
+++ b/src/mat/4C_mat_inelastic_defgrad_factors.hpp
@@ -654,6 +654,32 @@ namespace Mat
 
     virtual void unpack_inelastic(Core::Communication::UnpackBuffer& data) = 0;
 
+    /*!
+     * @brief Register names of the internal data that should be saved during runtime output
+     *
+     * @param[out] name_and_size Unordered map of names of the data with the respective vector size
+     */
+    virtual void register_output_data_names(
+        std::unordered_map<std::string, int>& names_and_size) const
+    {
+    }
+
+    /*!
+     * @brief Evaluate internal data for every Gauss point saved for output during runtime
+     * output
+     *
+     * @param[in] name  Name of the data to export
+     * @param[out] data NUMGPxNUMDATA Matrix holding the data
+     *
+     * @return true if data is set by the material, otherwise false
+     */
+    virtual bool evaluate_output_data(
+        const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+    {
+      return false;
+    }
+
+
    private:
     /// material parameters
     Core::Mat::PAR::Parameter* params_;
@@ -1326,6 +1352,12 @@ namespace Mat
 
     void unpack_inelastic(Core::Communication::UnpackBuffer& buffer) override;
 
+    void register_output_data_names(
+        std::unordered_map<std::string, int>& names_and_size) const override;
+
+    bool evaluate_output_data(
+        const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const override;
+
     /*! @brief Evaluate the current state variables based on a given right Cauchy-Green
      * deformation tensor, given inverse plastic deformation gradient and given equivalent
      * plastic strain
diff --git a/src/mat/4C_mat_inelastic_defgrad_factors_service.cpp b/src/mat/4C_mat_inelastic_defgrad_factors_service.cpp
index 17d4fb0182..31aad243a6 100644
--- a/src/mat/4C_mat_inelastic_defgrad_factors_service.cpp
+++ b/src/mat/4C_mat_inelastic_defgrad_factors_service.cpp
@@ -170,12 +170,17 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
   current_plastic_strain.resize(1, 0.0);  // value irrelevant at this point
   last_substep_plastic_strain.resize(1, 0.0);
 
+  // update last_ and current_ values of the equivalent stress
+  last_equiv_stress.resize(1, 0.0);
+  current_equiv_stress.resize(1, 0.0);  // value irrelevant at this point
+
   // default values of the right CG tensor: unit tensor
   last_rightCG.resize(1, id3x3);
   current_rightCG.resize(1, id3x3);  // value irrelevant at this point
 
   // default value for the current deformation gradient: zero tensor \f$ \boldsymbol{0} f$ (to make
   // sure that the inverse inelastic deformation gradient is evaluated in the first method call)
+  last_defgrad.resize(1, Core::LinAlg::Matrix<3, 3>{id3x3});
   current_defgrad.resize(1, Core::LinAlg::Matrix<3, 3>{Core::LinAlg::Initialization::zero});
 }
 
@@ -190,6 +195,7 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
       "you attempt to set it to {}",
       last_plastic_strain.size(), numgp);
 
+
   // default values of the inverse plastic deformation gradient for ALL Gauss Points
   last_plastic_defgrad_inverse.resize(numgp, last_plastic_defgrad_inverse[0]);
   current_plastic_defgrad_inverse.resize(numgp,
@@ -205,7 +211,12 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
   last_rightCG.resize(numgp, last_rightCG[0]);
   current_rightCG.resize(numgp, last_rightCG[0]);  // value irrelevant at this point
 
+  // default values of the equivalent stress for ALL Gauss Points
+  last_equiv_stress.resize(numgp, last_equiv_stress[0]);
+  current_equiv_stress.resize(numgp, current_equiv_stress[0]);  // value irrelevant at this point
+
   // default values of the deformation gradient
+  last_defgrad.resize(numgp, last_defgrad[0]);
   current_defgrad.resize(numgp, current_defgrad[0]);
 }
 
@@ -224,10 +235,12 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
 void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities::update()
 {
   // update history variables for the next time step
+  last_defgrad = current_defgrad;
   last_rightCG = current_rightCG;
   last_plastic_defgrad_inverse = current_plastic_defgrad_inverse;
   last_substep_plastic_defgrad_inverse = current_plastic_defgrad_inverse;
   last_plastic_strain = current_plastic_strain;
+  last_equiv_stress = current_equiv_stress;
   last_substep_plastic_strain = current_plastic_strain;
 }
 
@@ -237,9 +250,11 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
 void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities::pack(
     Core::Communication::PackBuffer& data) const
 {
+  add_to_pack(data, last_defgrad);
   add_to_pack(data, last_rightCG);
   add_to_pack(data, last_plastic_defgrad_inverse);
   add_to_pack(data, last_plastic_strain);
+  add_to_pack(data, last_equiv_stress);
   add_to_pack(data, last_substep_plastic_defgrad_inverse);
   add_to_pack(data, last_substep_plastic_strain);
 }
@@ -250,9 +265,11 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
     Core::Communication::UnpackBuffer& buffer)
 {
   // extract last values
+  extract_from_pack(buffer, last_defgrad);
   extract_from_pack(buffer, last_rightCG);
   extract_from_pack(buffer, last_plastic_defgrad_inverse);
   extract_from_pack(buffer, last_plastic_strain);
+  extract_from_pack(buffer, last_equiv_stress);
   extract_from_pack(buffer, last_substep_plastic_defgrad_inverse);
   extract_from_pack(buffer, last_substep_plastic_strain);
 
@@ -263,6 +280,8 @@ void Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::TimeStepQuantities:
       last_plastic_defgrad_inverse[0]);  // value irrelevant
   current_plastic_strain.resize(last_plastic_strain.size(),
       last_plastic_strain[0]);  // value irrelevant
+  current_equiv_stress.resize(last_equiv_stress.size(),
+      last_equiv_stress[0]);  // value irrelevant
 
   // set evaluated deformation gradient to 0, to make sure that the inverse inelastic deformation
   // gradient is evaluated fully after the restart
diff --git a/src/mat/4C_mat_inelastic_defgrad_factors_service.hpp b/src/mat/4C_mat_inelastic_defgrad_factors_service.hpp
index afb3a88b0b..493dff910f 100644
--- a/src/mat/4C_mat_inelastic_defgrad_factors_service.hpp
+++ b/src/mat/4C_mat_inelastic_defgrad_factors_service.hpp
@@ -62,8 +62,14 @@ namespace Mat
                           ///< analytical evaluation of the linearization
       failed_solution_analytic_linearization,  ///< solution of the linear system in the analytical
                                                ///< linearization failed
-      failed_matrix_log_evaluation,            ///< failed evaluation of the matrix logarithm or its
-                                               ///< derivative
+      failed_computation_flow_resistance,  ///< failed in the computation of the flow resistance via
+                                           ///< time integration of the hardening-rate equation
+                                           ///(e.g., when using the Anand law)
+      failed_computation_flow_resistance_derivs,  ///< failed in the computation of the flow
+                                                  ///< resistance derivatives (e.g., when using the
+                                                  ///< Anand law)
+      failed_matrix_log_evaluation,   ///< failed evaluation of the matrix logarithm or its
+                                      ///< derivative
       failed_matrix_exp_evaluation,   ///< failed evaluation of the matrix exponential or its
                                       ///< derivative
       failed_right_cg_interpolation,  ///< failed interpolation of the right Cauchy-Green tensor
@@ -133,6 +139,13 @@ namespace Mat
       //! (equivalent) plastic strain at the last time step (for all Gauss points)
       std::vector<double> last_plastic_strain;
 
+      //! equivalent stress at the previous time instant (for all Gauss points)
+      std::vector<double> last_equiv_stress;
+
+
+      //! last (reduced) deformation gradient (for all Gauss points)
+      std::vector<Core::LinAlg::Matrix<3, 3>> last_defgrad;
+
       //! temporary variable, for which we store the right Cauchy-Green deformation tensor at each
       //! evaluation (used in order to update last_rightCG_ once outer NR converges) (for all Gauss
       //! points)
@@ -148,6 +161,10 @@ namespace Mat
       //! current plastic strain (for all Gauss points)
       std::vector<double> current_plastic_strain;
 
+      //! current equivalent stress (for all Gauss points)
+      std::vector<double> current_equiv_stress;
+
+
       //! inverse plastic deformation gradient at the last computed time instant (after the last
       //! converged substep)
       std::vector<Core::LinAlg::Matrix<3, 3>> last_substep_plastic_defgrad_inverse;
@@ -422,7 +439,17 @@ namespace Mat
       StateQuantityDerivEvalType eval_type;
     };
 
+    //! struct containing specific derivatives of the scalar plastic strain rate used in
+    //! InelasticDefgradTransvIsotropElastViscoplast
+    struct PlasticStrainRateDerivs
+    {
+      //! derivative with respect to the equivalent stress
+      double deriv_equiv_stress;
+
 
+      //! derivative with respect to the plastic strain
+      double deriv_plastic_strain;
+    };
 
   }  // namespace InelasticDefgradTransvIsotropElastViscoplastUtils
 
diff --git a/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.cpp b/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.cpp
index 7f3696d6f8..ba09bd982a 100644
--- a/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.cpp
+++ b/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.cpp
@@ -1077,6 +1077,20 @@ void Mat::MultiplicativeSplitDefgradElastHyper::set_concentration_gp(const doubl
     inelastic_->fac_def_grad_in()[p].second->set_concentration_gp(concentration);
 }
 
+void Mat::MultiplicativeSplitDefgradElastHyper::register_output_data_names(
+    std::unordered_map<std::string, int>& names_and_size) const
+{
+  // register output for the inelastic defgrad factors
+  inelastic_->register_output_data_names(names_and_size);
+}
+
+bool Mat::MultiplicativeSplitDefgradElastHyper::evaluate_output_data(
+    const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+{
+  // evaluate GP output for the inelastic defgrad factors
+  return inelastic_->evaluate_output_data(name, data);
+}
+
 /*--------------------------------------------------------------------*
  *--------------------------------------------------------------------*/
 void Mat::InelasticFactorsHandler::assign_to_source(
@@ -1229,5 +1243,34 @@ Mat::MultiplicativeSplitDefgradElastHyper::evaluate_kinematic_quantities(
   return quantities;
 }
 
+void Mat::InelasticFactorsHandler::register_output_data_names(
+    std::unordered_map<std::string, int>& names_and_size) const
+{
+  // loop over all inelastic contributions
+  for (int i = 0; i < num_inelastic_def_grad(); ++i)
+  {
+    facdefgradin_[i].second->register_output_data_names(names_and_size);
+  }
+}
+
+bool Mat::InelasticFactorsHandler::evaluate_output_data(
+    const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+{
+  // declare the number of inelastic defgrad factors for which we use Gauss
+  // point data
+  int num_inel_fac_GP_data = 0;
+
+  // loop over all inelastic contributions
+  for (int i = 0; i < num_inelastic_def_grad(); ++i)
+  {
+    if (facdefgradin_[i].second->evaluate_output_data(name, data))
+    {
+      ++num_inel_fac_GP_data;
+    };
+  }
+
+  return (num_inel_fac_GP_data > 0);
+}
+
 
 FOUR_C_NAMESPACE_CLOSE
diff --git a/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.hpp b/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.hpp
index 82f664ff29..6191f0b676 100644
--- a/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.hpp
+++ b/src/mat/4C_mat_multiplicative_split_defgrad_elasthyper.hpp
@@ -135,6 +135,24 @@ namespace Mat
     /// Unpack
     void unpack_inelastic(Core::Communication::UnpackBuffer& buffer);
 
+    /*!
+     * @brief Register names of the internal data that should be saved during runtime output
+     *
+     * @param[out] name_and_size Unordered map of names of the data with the respective vector size
+     */
+    void register_output_data_names(std::unordered_map<std::string, int>& names_and_size) const;
+
+    /*!
+     * @brief Evaluate internal data for every Gauss point saved for output during runtime
+     * output
+     *
+     * @param[in] name  Name of the data to export
+     * @param[out] data NUMGPxNUMDATA Matrix holding the data
+     *
+     * @return true if data is set by the material, otherwise false
+     */
+    bool evaluate_output_data(const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const;
+
    private:
     /// vector that holds pairs of inelastic contribution and respective source
     std::vector<std::pair<PAR::InelasticSource, std::shared_ptr<Mat::InelasticDefgradFactors>>>
@@ -287,6 +305,12 @@ namespace Mat
 
     void update() override;
 
+    void register_output_data_names(
+        std::unordered_map<std::string, int>& names_and_size) const override;
+
+    bool evaluate_output_data(
+        const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const override;
+
     /*!
      * @brief Evaluate off-diagonal stiffness matrix (required for monolithic algorithms)
      *
diff --git a/src/mat/vplast/4C_mat_vplast_law.hpp b/src/mat/vplast/4C_mat_vplast_law.hpp
index 8dff418f4f..1c19ded260 100644
--- a/src/mat/vplast/4C_mat_vplast_law.hpp
+++ b/src/mat/vplast/4C_mat_vplast_law.hpp
@@ -94,6 +94,22 @@ namespace Mat
       virtual double evaluate_stress_ratio(
           const double equiv_stress, const double equiv_plastic_strain) = 0;
 
+
+      /*!
+       * @brief Evaluates the yield stress, or more generally the flow resistance (to accommodate
+       * viscoplastic laws with and without yield surfaces in a common framework)
+       *
+       *
+       * @param[in] equiv_stress equivalent stress
+       * @param[in] equiv_plastic_strain accumulated plastic strain
+       * @param[in/out] err_status error status of the evaluation
+       *
+       */
+      virtual double compute_flow_resistance(const double equiv_stress,
+          const double equiv_plastic_strain,
+          Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status) = 0;
+
+
       /*!
        * @brief Evaluate the equivalent plastic strain rate \f$ \dot{\varepsilon}^{\text{p}} \f$ for
        * a given equivalent stress \f$ \overflow{\sigma} \f$ and a given plastic strain \f$
@@ -142,8 +158,9 @@ namespace Mat
        * @return Derivatives of the equivalent plastic strain rate w.r.t. the equivalent stress
        *         (element 0 of matrix) and the plastic strain (element 1 of matrix)
        */
-      virtual Core::LinAlg::Matrix<2, 1> evaluate_derivatives_of_plastic_strain_rate(
-          const double equiv_stress, const double equiv_plastic_strain, const double dt,
+      virtual InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs
+      evaluate_derivatives_of_plastic_strain_rate(const double equiv_stress,
+          const double equiv_plastic_strain, const double dt,
           const double max_plastic_strain_deriv_incr,
           Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status,
           const bool update_hist_var = true) = 0;
@@ -164,9 +181,10 @@ namespace Mat
        * @brief Pre-evaluation, intended to be used for stuff that has to be done only once per
        *        evaluate()
        *
+       * @param[in] params  parameter list
        * @param[in] gp      Current Gauss point
        */
-      virtual void pre_evaluate(int gp) = 0;
+      virtual void pre_evaluate(const Teuchos::ParameterList& params, int gp) { gp_ = gp; };
 
       /*!
        * @brief Update history variables of the viscoplasticity law for next time step
@@ -185,6 +203,38 @@ namespace Mat
 
       virtual void unpack_viscoplastic_law(Core::Communication::UnpackBuffer& buffer) = 0;
 
+
+      /*!
+       * @brief Register names of the internal data that should be saved during runtime output
+       *
+       * @param[out] name_and_size Unordered map of names of the data with the respective vector
+       * size
+       */
+      virtual void register_output_data_names(
+          std::unordered_map<std::string, int>& names_and_size) const
+      {
+      }
+
+      /*!
+       * @brief Evaluate internal data for every Gauss point saved for output during runtime
+       * output
+       *
+       * @param[in] name  Name of the data to export
+       * @param[out] data NUMGPxNUMDATA Matrix holding the data
+       *
+       * @return true if data is set by the material, otherwise false
+       */
+      virtual bool evaluate_output_data(
+          const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+      {
+        return false;
+      }
+
+     protected:
+      /// Gauss point index
+      int gp_;
+
+
      private:
       /// material parameters
       Core::Mat::PAR::Parameter* params_;
diff --git a/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.cpp b/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.cpp
index 6842b9a461..9f70952631 100644
--- a/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.cpp
+++ b/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.cpp
@@ -12,10 +12,12 @@
 #include "4C_mat_inelastic_defgrad_factors_service.hpp"
 #include "4C_mat_par_bundle.hpp"
 #include "4C_mat_vplast_law.hpp"
+#include "4C_utils_exceptions.hpp"
 
 #include <cmath>
 #include <utility>
 
+
 FOUR_C_NAMESPACE_OPEN
 
 
@@ -30,8 +32,18 @@ Mat::Viscoplastic::PAR::ReformulatedJohnsonCook::ReformulatedJohnsonCook(
       strain_rate_exp_fac_(matdata.parameters.get<double>("STRAIN_RATE_EXP_FAC")),
       init_yield_strength_(matdata.parameters.get<double>("INIT_YIELD_STRENGTH")),
       isotrop_harden_prefac_(matdata.parameters.get<double>("ISOTROP_HARDEN_PREFAC")),
-      isotrop_harden_exp_(matdata.parameters.get<double>("ISOTROP_HARDEN_EXP"))
+      isotrop_harden_exp_(matdata.parameters.get<double>("ISOTROP_HARDEN_EXP")),
+      ref_temperature_(matdata.parameters.get<double>("REF_TEMPERATURE")),
+      melt_temperature_(matdata.parameters.get<double>("MELT_TEMPERATURE")),
+      temperature_sens_(matdata.parameters.get<double>("TEMPERATURE_SENS"))
 {
+  // consistency checks for the temperature
+  FOUR_C_ASSERT_ALWAYS(ref_temperature_ > 0.0,
+      "Reference temperature must be > 0.0: current value: {}", ref_temperature_);
+  FOUR_C_ASSERT_ALWAYS(melt_temperature_ > ref_temperature_,
+      "Melting temperature must be > reference temperature: current values: {} (T_m) and {} "
+      "(T_ref)",
+      melt_temperature_, ref_temperature_);
 }
 
 /*--------------------------------------------------------------------*
@@ -45,20 +57,77 @@ Mat::Viscoplastic::ReformulatedJohnsonCook::ReformulatedJohnsonCook(
 {
 }
 
+
+/*--------------------------------------------------------------------*
+ *--------------------------------------------------------------------*/
+void Mat::Viscoplastic::ReformulatedJohnsonCook::pre_evaluate(
+    const Teuchos::ParameterList& params, int gp)
+{
+  // call pre_evaluate of base class
+  Mat::Viscoplastic::Law::pre_evaluate(params, gp);
+
+
+  // get temperature factors
+  double T = parameter()->ref_temperature();
+  if (params.isParameter("temperature"))
+  {
+    T = params.get<double>("temperature");
+  }
+  const double T_ref = parameter()->ref_temperature();
+  const double T_melt = parameter()->melt_temperature();
+  const double M = parameter()->temperature_sens();
+
+  // consistency checks for the temperature
+  FOUR_C_ASSERT_ALWAYS(T > 0.0 && T < T_melt,
+      "Temperature must be > 0.0 and < melting temperature: current value: {} (T) and {} (T_m)", T,
+      T_melt);
+
+  // set temperature ratio
+  temperature_ratio_ = 1.0;
+  if (T != T_ref)
+  {
+    FOUR_C_ASSERT_ALWAYS(T_ref < T_melt,
+        "You specified the reference temperature {} >= melting temperature: {}! This cannot be "
+        "currently resolved by the Reformulated Johnson-Cook law!",
+        T_ref, T_melt);
+    temperature_ratio_ =
+        1.0 - (std::pow(T, M) - std::pow(T_ref, M)) / (std::pow(T_melt, M) - std::pow(T_ref, M));
+  }
+}
+
 /*--------------------------------------------------------------------*
  *--------------------------------------------------------------------*/
 double Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_stress_ratio(
     const double equiv_stress, const double equiv_plastic_strain)
 {
   // extract yield strength from the plastic strain and the material parameters
+  InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType yield_strength_err_status{
+      InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors};
   const double yield_strength =
-      (parameter()->init_yield_strength() +
-          parameter()->isotrop_harden_prefac() *
-              std::pow(equiv_plastic_strain, parameter()->isotrop_harden_exp()));
+      compute_flow_resistance(equiv_stress, equiv_plastic_strain, yield_strength_err_status);
+  FOUR_C_ASSERT_ALWAYS(yield_strength_err_status ==
+                           InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors,
+      "Something went wrong when evaluating the yield stress: error = {}",
+      EnumTools::enum_name(yield_strength_err_status));
 
   return equiv_stress / yield_strength;
 }
 
+
+/*--------------------------------------------------------------------*
+ *--------------------------------------------------------------------*/
+double Mat::Viscoplastic::ReformulatedJohnsonCook::compute_flow_resistance(
+    const double equiv_stress, const double equiv_plastic_strain,
+    Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status)
+{
+  // extract yield strength from the plastic strain and the material parameters
+  return (parameter()->init_yield_strength() * temperature_ratio_ +
+          parameter()->isotrop_harden_prefac() * temperature_ratio_ *
+              std::pow(equiv_plastic_strain, parameter()->isotrop_harden_exp()));
+}
+
+
+
 /*--------------------------------------------------------------------*
  *--------------------------------------------------------------------*/
 double Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_plastic_strain_rate(
@@ -80,6 +149,32 @@ double Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_plastic_strain_rate(
   // compute the viscoplastic strain rate; first we set it to 0
   double equiv_plastic_strain_rate = 0.0;
 
+  // save yield strength
+  if (update_hist_var)
+  {
+    // GP index safeguard
+    FOUR_C_ASSERT_ALWAYS(
+        static_cast<size_t>(gp_) < time_step_quantities_.current_yield_strength_.size(),
+        "The current gp index is {} while the stored number of GP within Reformulated Johnson - "
+        "Cook is {} ",
+        gp_, time_step_quantities_.current_yield_strength_.size());
+    InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType yield_strength_err_status =
+        InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors;
+    time_step_quantities_.current_yield_strength_[gp_] =
+        compute_flow_resistance(equiv_stress, equiv_plastic_strain, yield_strength_err_status);
+    FOUR_C_ASSERT_ALWAYS(
+        yield_strength_err_status ==
+            InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors,
+        "Computing yield strength during the evaluation of the plastic strain rate has failed with "
+        "error {}",
+        EnumTools::enum_name(yield_strength_err_status));
+  }
+
+  // verify whether the maximum plastic strain increment is larger than 0, since we will take its
+  // logarithm
+  FOUR_C_ASSERT_ALWAYS(max_plastic_strain_incr > 0.0,
+      "Maximum plastic strain increment must be > 0: current value = {}", max_plastic_strain_incr);
+
   // stress ratio
   double stress_ratio = evaluate_stress_ratio(equiv_stress, equiv_plastic_strain);
 
@@ -106,18 +201,21 @@ double Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_plastic_strain_rate(
 
 /*--------------------------------------------------------------------*
  *--------------------------------------------------------------------*/
-Core::LinAlg::Matrix<2, 1>
+Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs
 Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_derivatives_of_plastic_strain_rate(
     const double equiv_stress, const double equiv_plastic_strain, const double dt,
     const double max_plastic_strain_deriv_incr, ViscoplastErrorType& err_status,
     const bool update_hist_var)
 {
+  // declare derivatives to be computed
+  double deriv_equiv_stress{0.0};
+  double deriv_plastic_strain{0.0};
+
   // first set error status to "no errors"
   err_status = ViscoplastErrorType::no_errors;
 
   // used equivalent plastic strain
   double used_equiv_plastic_strain = equiv_plastic_strain;
-
   // check whether the plastic strain is less than a set value (singularity in the derivatives
   // below)
   if (std::abs(equiv_plastic_strain) < 1.0e-16)
@@ -130,12 +228,20 @@ Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_derivatives_of_plastic_stra
   if (equiv_plastic_strain < 0.0)
   {
     err_status = ViscoplastErrorType::negative_plastic_strain;
-    return Core::LinAlg::Matrix<2, 1>{Core::LinAlg::Initialization::zero};
+    return InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs{
+        .deriv_equiv_stress = 0.0, .deriv_plastic_strain = 0.0};
   }
 
   // extraction of the yield strength from the plastic strain and the material parameters
+  InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType yield_strength_err_status{
+      InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors};
   const double yield_strength =
-      const_pars_.sigma_Y0 + const_pars_.B * std::pow(used_equiv_plastic_strain, const_pars_.N);
+      compute_flow_resistance(equiv_stress, used_equiv_plastic_strain, yield_strength_err_status);
+  FOUR_C_ASSERT_ALWAYS(yield_strength_err_status ==
+                           InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors,
+      "Something went wrong in the computation of the yield strength in the plastic strain rate "
+      "derivative evaluation: error {}",
+      EnumTools::enum_name(yield_strength_err_status));
   const double log_yield_strength = std::log(yield_strength);
   const double inv_yield_strength = 1.0 / yield_strength;
 
@@ -149,41 +255,114 @@ Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_derivatives_of_plastic_stra
 
   // computation of derivatives
 
-  // first we set derivatives to 0
-  Core::LinAlg::Matrix<2, 1> equiv_plastic_strain_rate_ders(Core::LinAlg::Initialization::zero);
-
   // then we check the yield condition
-  if (evaluate_stress_ratio(equiv_stress, used_equiv_plastic_strain) >= 1.0)
+  if (evaluate_stress_ratio(equiv_stress, equiv_plastic_strain) >= 1.0)
   {
     // compute first the logarithms of our derivatives (try to avoid overflow!)
     double log_deriv_sigma = const_pars_.log_p_e +
                              const_pars_.e * (equiv_stress * inv_yield_strength - 1.0) -
                              log_yield_strength;
-    double log_deriv_eps = const_pars_.log_p_e +
-                           const_pars_.e * (equiv_stress * inv_yield_strength - 1.0) +
-                           log_equiv_stress - 2.0 * log_yield_strength + const_pars_.log_B_N +
-                           (const_pars_.N - 1.0) * log_equiv_plastic_strain;
 
-    // check overflow error using these logarithms
-    if (max_plastic_strain_deriv_incr <= 0.0)
+    // verify whether the maximum plastic strain derivative increment is larger than 0, since we
+    // will take its logarithm
+    FOUR_C_ASSERT_ALWAYS(max_plastic_strain_deriv_incr > 0.0,
+        "Maximum plastic strain derivative increment must be > 0: current value = {}",
+        max_plastic_strain_deriv_incr);
+
+    // perfect plasticity
+    if (const_pars_.is_perfect_plasticity)
     {
-      err_status = ViscoplastErrorType::overflow_error;
-      return Core::LinAlg::Matrix<2, 1>{Core::LinAlg::Initialization::zero};
+      // check overflow error using these logarithms
+      double log_max_plastic_strain_deriv_value = std::log(max_plastic_strain_deriv_incr);
+      if ((log_dt + log_deriv_sigma > log_max_plastic_strain_deriv_value))
+      {
+        err_status = ViscoplastErrorType::failed_computation_flow_resistance_derivs;
+        return InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs{
+            .deriv_equiv_stress = 0.0, .deriv_plastic_strain = 0.0};
+      }
+
+      // compute the exact derivatives using these logarithms
+      deriv_equiv_stress = std::exp(log_deriv_sigma);
+      deriv_plastic_strain = 0.0;
     }
-    double log_max_plastic_strain_deriv_value = std::log(max_plastic_strain_deriv_incr);
-    if ((log_dt + log_deriv_sigma > log_max_plastic_strain_deriv_value) &&
-        (log_dt + log_deriv_eps > log_max_plastic_strain_deriv_value))
+    // hardening case
+    else
     {
-      err_status = ViscoplastErrorType::overflow_error;
-      return Core::LinAlg::Matrix<2, 1>{Core::LinAlg::Initialization::zero};
+      const double log_deriv_eps =
+          const_pars_.log_p_e + const_pars_.e * (equiv_stress * inv_yield_strength - 1.0) +
+          log_equiv_stress - 2.0 * log_yield_strength + const_pars_.log_B_N +
+          (const_pars_.N - 1.0) * log_equiv_plastic_strain;
+      // check overflow error using these logarithms
+      double log_max_plastic_strain_deriv_value = std::log(max_plastic_strain_deriv_incr);
+      if ((log_dt + log_deriv_sigma > log_max_plastic_strain_deriv_value) &&
+          (log_dt + log_deriv_eps > log_max_plastic_strain_deriv_value))
+      {
+        err_status = ViscoplastErrorType::failed_computation_flow_resistance_derivs;
+        return InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs{
+            .deriv_equiv_stress = 0.0, .deriv_plastic_strain = 0.0};
+      }
+
+      // compute the exact derivatives using these logarithms
+      deriv_equiv_stress = std::exp(log_deriv_sigma);
+      deriv_plastic_strain = -std::exp(log_deriv_eps);
+    }
+  }
+
+  return InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs{
+      .deriv_equiv_stress = deriv_equiv_stress, .deriv_plastic_strain = deriv_plastic_strain};
+}
+
+void Mat::Viscoplastic::ReformulatedJohnsonCook::setup(const int numgp,
+    const Discret::Elements::Fibers& fibers,
+    const std::optional<Discret::Elements::CoordinateSystem>& coord_system)
+{
+  time_step_quantities_.current_yield_strength_.resize(numgp, parameter()->init_yield_strength());
+}
+
+void Mat::Viscoplastic::ReformulatedJohnsonCook::register_output_data_names(
+    std::unordered_map<std::string, int>& names_and_size) const
+{
+  names_and_size["yield_strength"] = 1;
+}
+
+bool Mat::Viscoplastic::ReformulatedJohnsonCook::evaluate_output_data(
+    const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const
+{
+  if (name == "yield_strength")
+  {
+    for (int gp = 0; gp < static_cast<int>(time_step_quantities_.current_yield_strength_.size());
+        ++gp)
+    {
+      data(gp, 0) = time_step_quantities_.current_yield_strength_[gp];
     }
+    return true;
+  }
+
+  return false;
+}
 
-    // compute the exact derivatives using these logarithms
-    equiv_plastic_strain_rate_ders(0, 0) = std::exp(log_deriv_sigma);
-    equiv_plastic_strain_rate_ders(1, 0) = -std::exp(log_deriv_eps);
+void Mat::Viscoplastic::ReformulatedJohnsonCook::pack_viscoplastic_law(
+    Core::Communication::PackBuffer& data) const
+{
+  // pack relevant variables
+  if (parameter() != nullptr)
+  {
+    // we need to pack this current value since it also saves the number
+    // of Gauss points
+    add_to_pack(data, time_step_quantities_.current_yield_strength_);
   }
+}
 
-  return equiv_plastic_strain_rate_ders;
+void Mat::Viscoplastic::ReformulatedJohnsonCook::unpack_viscoplastic_law(
+    Core::Communication::UnpackBuffer& buffer)
+{
+  // NOTE: factory method is called during assign_to_source in the unpack method of the
+  // multiplicative split framework --> the param class of the viscoplastic law is created, we only
+  // need to unpack the history variables
+  if (parameter() != nullptr)
+  {
+    extract_from_pack(buffer, time_step_quantities_.current_yield_strength_);
+  }
 }
 
 FOUR_C_NAMESPACE_CLOSE
diff --git a/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.hpp b/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.hpp
index 04bb5a147c..8e18d3fd15 100644
--- a/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.hpp
+++ b/src/mat/vplast/4C_mat_vplast_reform_johnsoncook.hpp
@@ -58,6 +58,12 @@ namespace Mat
         [[nodiscard]] double isotrop_harden_prefac() const { return isotrop_harden_prefac_; };
         //! get exponent of the isotropic hardening stress \f$ n \f$
         [[nodiscard]] double isotrop_harden_exp() const { return isotrop_harden_exp_; };
+        //! get reference temperature
+        [[nodiscard]] double ref_temperature() const { return ref_temperature_; };
+        //! get melting temperature
+        [[nodiscard]] double melt_temperature() const { return melt_temperature_; };
+        //! get temperature sensitivity factor
+        [[nodiscard]] double temperature_sens() const { return temperature_sens_; };
 
        private:
         //! strain rate prefactor \f$ \dot{P}_0 \f$
@@ -75,6 +81,15 @@ namespace Mat
 
         //! exponent of the isotropic hardening stress \f$ n \f$
         const double isotrop_harden_exp_;
+
+        //! reference temperature \f$ T_{\mathrm{ref}} \f$
+        const double ref_temperature_;
+
+        //! melting temperature \f$ T_{\mathrm{melt}} \f$
+        const double melt_temperature_;
+
+        //! temperature sensitivity \f$ M \f$
+        const double temperature_sens_;
       };
     }  // namespace PAR
 
@@ -104,35 +119,48 @@ namespace Mat
       double evaluate_stress_ratio(
           const double equiv_stress, const double equiv_plastic_strain) override;
 
+      double compute_flow_resistance(const double equiv_stress, const double equiv_plastic_strain,
+          Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status) override;
 
       double evaluate_plastic_strain_rate(const double equiv_stress,
           const double equiv_plastic_strain, const double dt, const double max_plastic_strain_incr,
           Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status,
-          const bool update_hist_var = true) override;
+          const bool update_hist_var) override;
 
-      Core::LinAlg::Matrix<2, 1> evaluate_derivatives_of_plastic_strain_rate(
-          const double equiv_stress, const double equiv_plastic_strain, const double dt,
+      InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs
+      evaluate_derivatives_of_plastic_strain_rate(const double equiv_stress,
+          const double equiv_plastic_strain, const double dt,
           const double max_plastic_strain_deriv_incr,
           Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType& err_status,
-          const bool update_hist_var = true) override;
+          const bool update_hist_var) override;
 
       void setup(const int numgp, const Discret::Elements::Fibers& fibers,
-          const std::optional<Discret::Elements::CoordinateSystem>& coord_system) override {};
+          const std::optional<Discret::Elements::CoordinateSystem>& coord_system) override;
+
 
-      void pre_evaluate(int gp) override {};
+      void pre_evaluate(const Teuchos::ParameterList& params, int gp) override;
 
       void update() override {};
 
       void update_gp_state(int gp) override {};
 
-      void pack_viscoplastic_law(Core::Communication::PackBuffer& data) const override {};
+      void pack_viscoplastic_law(Core::Communication::PackBuffer& data) const override;
 
-      void unpack_viscoplastic_law(Core::Communication::UnpackBuffer& buffer) override {};
+      void unpack_viscoplastic_law(Core::Communication::UnpackBuffer& buffer) override;
+
+      void register_output_data_names(
+          std::unordered_map<std::string, int>& names_and_size) const override;
+
+      bool evaluate_output_data(
+          const std::string& name, Core::LinAlg::SerialDenseMatrix& data) const override;
 
      private:
       /// struct containing constant parameters to be evaluated only once
       struct ConstPars
       {
+        /// perfect plasticity? \f$ B * N == 0 \f$
+        bool is_perfect_plasticity;
+
         /// prefactor \f$ P \f$ of the plastic strain rate
         double p;
 
@@ -160,20 +188,56 @@ namespace Mat
         /// constructor
         ConstPars(const double prefac, const double expon, const double harden_prefac,
             const double harden_expon, const double initial_yield_strength)
-            : p(prefac),
+            : is_perfect_plasticity(harden_prefac == 0.0 || harden_expon == 0.0),
+              p(prefac),
               log_p(std::log(p)),
               e(expon),
               log_p_e(std::log(prefac * expon)),
               B(harden_prefac),
               N(harden_expon),
-              log_B_N(std::log(harden_prefac * harden_expon)),
+              log_B_N(is_perfect_plasticity ? 0.0 : set_log_b_n(harden_prefac, harden_expon)),
               sigma_Y0(initial_yield_strength)
         {
+          FOUR_C_ASSERT_ALWAYS(B >= 0.0,
+              "The current Reformulated Johnson-Cook viscoplastic law only enables hardening "
+              "prefactors >= 0: current specification: {}",
+              B);
+          FOUR_C_ASSERT_ALWAYS(N >= 0.0,
+              "The current Reformulated Johnson-Cook viscoplastic law only enables hardening "
+              "exponents >= 0: current specification: {}",
+              N);
+        }
+
+        /*!
+         * @brief Computes log(B * N) for the Johnson-Cook hardening law, assuming B * N > 0.0
+         *
+         *
+         *@param[in] B Hardening prefactor.
+         *@param[in] Hardening exponent.
+         *@return logarithm log(B * N)
+         */
+        double set_log_b_n(const double B, const double N)
+        {
+          FOUR_C_ASSERT_ALWAYS(
+              B * N > 0.0, "You should not be here when setting B * N: {} * {} <= 0.0", B, N);
+          return std::log(B * N);
         }
       };
 
       /// instance of ConstPars struct
       const ConstPars const_pars_;
+
+      /// temperature ratio \f$ D_T = 1 - \frac{T^M - T_{\mathrm{ref}}^M}{T_{\mathrm{melt}}^M -
+      /// T_{\mathrm{ref}}^M} \f$
+      double temperature_ratio_ = 1.0;
+
+      //! struct containing quantities at different time points (i.e., "current" at \f[ t_{n+1} \f])
+      struct TimeStepQuantities
+      {
+        //! yield strength (for all Gauss points)
+        std::vector<double> current_yield_strength_;
+      };
+      TimeStepQuantities time_step_quantities_;
     };
 
   }  // namespace Viscoplastic
diff --git a/tests/input_files/mat_iso_viscoplast_refJC_log_substep.4C.yaml b/tests/input_files/mat_iso_viscoplast_refJC_log_substep.4C.yaml
index 6b54bec580..4e6b283fcc 100644
--- a/tests/input_files/mat_iso_viscoplast_refJC_log_substep.4C.yaml
+++ b/tests/input_files/mat_iso_viscoplast_refJC_log_substep.4C.yaml
@@ -59,6 +59,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 5
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/tests/input_files/mat_transviso_viscoplast_refJC_log_substep.4C.yaml b/tests/input_files/mat_transviso_viscoplast_refJC_log_substep.4C.yaml
index 2fbfd1b42e..ef16e1d02c 100644
--- a/tests/input_files/mat_transviso_viscoplast_refJC_log_substep.4C.yaml
+++ b/tests/input_files/mat_transviso_viscoplast_refJC_log_substep.4C.yaml
@@ -62,6 +62,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 5
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/tests/input_files/mat_transviso_viscoplast_refJC_standard_substep.4C.yaml b/tests/input_files/mat_transviso_viscoplast_refJC_standard_substep.4C.yaml
index 63ab0466df..f74f4e4943 100644
--- a/tests/input_files/mat_transviso_viscoplast_refJC_standard_substep.4C.yaml
+++ b/tests/input_files/mat_transviso_viscoplast_refJC_standard_substep.4C.yaml
@@ -62,6 +62,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 5
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/tests/input_files/poro_new_solid_nitsche_contact3D_new_struc_vplast.4C.yaml b/tests/input_files/poro_new_solid_nitsche_contact3D_new_struc_vplast.4C.yaml
index b74d21959a..23acc66344 100644
--- a/tests/input_files/poro_new_solid_nitsche_contact3D_new_struc_vplast.4C.yaml
+++ b/tests/input_files/poro_new_solid_nitsche_contact3D_new_struc_vplast.4C.yaml
@@ -96,6 +96,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 24
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/tests/input_files/poro_new_solid_nitsche_contact3D_twoelem_new_struc_vplast.4C.yaml b/tests/input_files/poro_new_solid_nitsche_contact3D_twoelem_new_struc_vplast.4C.yaml
index d17bbc858f..5231b286eb 100644
--- a/tests/input_files/poro_new_solid_nitsche_contact3D_twoelem_new_struc_vplast.4C.yaml
+++ b/tests/input_files/poro_new_solid_nitsche_contact3D_twoelem_new_struc_vplast.4C.yaml
@@ -95,6 +95,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 24
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/tests/input_files/poro_solid_nitsche_contact3D_twoelem_twosided_new_struc_vplast.4C.yaml b/tests/input_files/poro_solid_nitsche_contact3D_twoelem_twosided_new_struc_vplast.4C.yaml
index b43780bac9..2a9699aa15 100644
--- a/tests/input_files/poro_solid_nitsche_contact3D_twoelem_twosided_new_struc_vplast.4C.yaml
+++ b/tests/input_files/poro_solid_nitsche_contact3D_twoelem_twosided_new_struc_vplast.4C.yaml
@@ -88,6 +88,9 @@ MATERIALS:
       INIT_YIELD_STRENGTH: 792
       ISOTROP_HARDEN_PREFAC: 510
       ISOTROP_HARDEN_EXP: 0.26
+      REF_TEMPERATURE: 293.0
+      TEMPERATURE_SENS: 1.03
+      MELT_TEMPERATURE: 1000.0
   - MAT: 24
     ELAST_CoupTransverselyIsotropic:
       ALPHA: 1
diff --git a/unittests/mat/4C_inelastic_defgrad_factors_test.cpp b/unittests/mat/4C_inelastic_defgrad_factors_test.cpp
index 154764e24f..527bcb928e 100644
--- a/unittests/mat/4C_inelastic_defgrad_factors_test.cpp
+++ b/unittests/mat/4C_inelastic_defgrad_factors_test.cpp
@@ -385,6 +385,12 @@ namespace
       viscoplastic_law_reformulated_Johnson_Cook_data.add("INIT_YIELD_STRENGTH", 2.0e4);
       viscoplastic_law_reformulated_Johnson_Cook_data.add("ISOTROP_HARDEN_PREFAC", 5.0e3);
       viscoplastic_law_reformulated_Johnson_Cook_data.add("ISOTROP_HARDEN_EXP", 0.2);
+      viscoplastic_law_reformulated_Johnson_Cook_data.add("REF_TEMPERATURE", 293.0);
+      viscoplastic_law_reformulated_Johnson_Cook_data.add("MELT_TEMPERATURE", 1793.0);
+      viscoplastic_law_reformulated_Johnson_Cook_data.add("TEMPERATURE_SENS", 1.03);
+
+
+
       // add material to problem instance
       problem.materials()->insert(400,
           Mat::make_parameter(400, Core::Materials::MaterialType::mvl_reformulated_Johnson_Cook,
diff --git a/unittests/mat/vplast/4C_vplast_reform_johnsoncook_test.cpp b/unittests/mat/vplast/4C_vplast_reform_johnsoncook_test.cpp
index 91bac63c36..de037f5084 100644
--- a/unittests/mat/vplast/4C_vplast_reform_johnsoncook_test.cpp
+++ b/unittests/mat/vplast/4C_vplast_reform_johnsoncook_test.cpp
@@ -10,6 +10,7 @@
 #include "4C_global_data.hpp"
 #include "4C_linalg_fixedsizematrix.hpp"
 #include "4C_mat_inelastic_defgrad_factors.hpp"
+#include "4C_mat_inelastic_defgrad_factors_service.hpp"
 #include "4C_mat_material_factory.hpp"
 #include "4C_mat_vplast_law.hpp"
 #include "4C_mat_vplast_reform_johnsoncook.hpp"
@@ -38,6 +39,9 @@ namespace
       vplast_law_reformulated_JC_data.add("INIT_YIELD_STRENGTH", 792.0);
       vplast_law_reformulated_JC_data.add("ISOTROP_HARDEN_PREFAC", 510.0);
       vplast_law_reformulated_JC_data.add("ISOTROP_HARDEN_EXP", 0.26);
+      vplast_law_reformulated_JC_data.add("REF_TEMPERATURE", 293.0);
+      vplast_law_reformulated_JC_data.add("MELT_TEMPERATURE", 1793.0);
+      vplast_law_reformulated_JC_data.add("TEMPERATURE_SENS", 1.03);
       params_vplast_law_reformulated_JC_ =
           std::dynamic_pointer_cast<Mat::Viscoplastic::PAR::ReformulatedJohnsonCook>(
               std::shared_ptr(Mat::make_parameter(1,
@@ -50,8 +54,13 @@ namespace
       int numgp = 8;  // HEX8 element, although not really relevant for the tested methods
       vplast_law_reformulated_JC_->setup(numgp, {}, {});
 
+      // parameter list
+      Teuchos::ParameterList param_list{};
+      param_list.set<double>("temperature", 293);
+
+
       // call pre_evaluate
-      vplast_law_reformulated_JC_->pre_evaluate(0);
+      vplast_law_reformulated_JC_->pre_evaluate(param_list, 0);
     }
 
     // equivalent stress
@@ -64,7 +73,8 @@ namespace
     double plastic_strain_rate_reformulated_JC_solution_;
     // reference solution for the plastic strain rate derivatives, w.r.t. equivalent stress and
     // plastic strain (ReformulatedJohnsonCook)
-    Core::LinAlg::Matrix<2, 1> deriv_plastic_strain_rate_reformulated_JC_solution_;
+    Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs
+        deriv_plastic_strain_rate_reformulated_JC_solution_;
     // pointer to ReformulatedJohnsonCook
     std::shared_ptr<Mat::Viscoplastic::ReformulatedJohnsonCook> vplast_law_reformulated_JC_;
     // pointer to parameters of ReformulatedJohnsonCook
@@ -113,8 +123,8 @@ namespace
   TEST_F(ReformJohnsonCookTest, TestEvaluatePlasticStrainRateDerivatives)
   {
     // set reference solution
-    deriv_plastic_strain_rate_reformulated_JC_solution_(0, 0) = 1889.49890189991;
-    deriv_plastic_strain_rate_reformulated_JC_solution_(1, 0) = -47431778.9968811;
+    deriv_plastic_strain_rate_reformulated_JC_solution_.deriv_equiv_stress = 1889.49890189991;
+    deriv_plastic_strain_rate_reformulated_JC_solution_.deriv_plastic_strain = -47431778.9968811;
 
 
     // declare error status
@@ -122,16 +132,19 @@ namespace
         Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors;
 
     // compute solution from the viscoplasticity law
-    Core::LinAlg::Matrix<2, 1> deriv_plastic_strain_rate_reformulated_JC =
-        vplast_law_reformulated_JC_->evaluate_derivatives_of_plastic_strain_rate(
-            equiv_stress_, equiv_plastic_strain_, 1.0, std::exp(30.0), err_status, false);
+    Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::PlasticStrainRateDerivs
+        deriv_plastic_strain_rate_reformulated_JC =
+            vplast_law_reformulated_JC_->evaluate_derivatives_of_plastic_strain_rate(
+                equiv_stress_, equiv_plastic_strain_, 1.0, std::exp(30.0), err_status, false);
 
     if (err_status != Mat::InelasticDefgradTransvIsotropElastViscoplastUtils::ErrorType::no_errors)
       FOUR_C_THROW("Error encountered during testing of TestEvaluatePlasticStrainRateDerivatives");
 
     // compare solutions
-    FOUR_C_EXPECT_NEAR(deriv_plastic_strain_rate_reformulated_JC_solution_,
-        deriv_plastic_strain_rate_reformulated_JC, 1.0e-6);
+    EXPECT_NEAR(deriv_plastic_strain_rate_reformulated_JC_solution_.deriv_equiv_stress,
+        deriv_plastic_strain_rate_reformulated_JC.deriv_equiv_stress, 1.0e-6);
+    EXPECT_NEAR(deriv_plastic_strain_rate_reformulated_JC_solution_.deriv_plastic_strain,
+        deriv_plastic_strain_rate_reformulated_JC.deriv_plastic_strain, 1.0e-6);
   }
 
 }  // namespace