diff --git a/src/JacRes.cpp b/src/JacRes.cpp
index c10d6202..7357125e 100644
--- a/src/JacRes.cpp
+++ b/src/JacRes.cpp
@@ -1930,9 +1930,11 @@ PetscErrorCode JacResCopyPres(JacRes *jr, Vec x)
 	PetscFunctionReturn(0);
 }
 //---------------------------------------------------------------------------
-PetscErrorCode JacResInitPres(JacRes *jr)
+PetscErrorCode JacResInitPres(JacRes *jr,TSSol *ts)
+
 {
 	FDSTAG            *fs;
+	
 	BCCtx             *bc;
 	SolVarCell        *svCell;
 	const PetscScalar *p;
@@ -1949,7 +1951,7 @@ PetscErrorCode JacResInitPres(JacRes *jr)
 	fixPhase = bc->fixPhase;
 
 	// check activation
-	if(!bc->initPres) PetscFunctionReturn(0);
+	if(!bc->initPres || ts->istep>0) PetscFunctionReturn(0);
 
 	// get grid coordinate bounds in z-direction
 	ierr = FDSTAGGetGlobalBox(fs, NULL, NULL, &bz, NULL, NULL, &ez); CHKERRQ(ierr);
@@ -1999,7 +2001,7 @@ PetscErrorCode JacResInitPres(JacRes *jr)
 	PetscFunctionReturn(0);
 }
 //---------------------------------------------------------------------------
-PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx)
+PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx,TSSol *ts)
 {
 	FDSTAG            *fs;
 	SolVarCell        *svCell;
@@ -2015,7 +2017,7 @@ PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx)
 	PetscFunctionBeginUser;
 
 	// check activation
-	if(!jr->ctrl.initLithPres) PetscFunctionReturn(0);
+	if(!jr->ctrl.initLithPres || ts->istep >0) PetscFunctionReturn(0);
 
 	// print
 	PrintStart(&t, "Initializing pressure with lithostatic pressure", NULL);
diff --git a/src/JacRes.h b/src/JacRes.h
index 994b3f76..833c0113 100644
--- a/src/JacRes.h
+++ b/src/JacRes.h
@@ -313,10 +313,10 @@ PetscErrorCode JacResCopyVel(JacRes *jr, Vec x);
 PetscErrorCode JacResCopyPres(JacRes *jr, Vec x);
 
 // initialize pressure
-PetscErrorCode JacResInitPres(JacRes *jr);
+PetscErrorCode JacResInitPres(JacRes *jr,TSSol *ts);
 
 // initialize pressure to lithostatic pressure
-PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx);
+PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx, TSSol *ts);
 
 // copy residuals from local to global vectors, enforce boundary constraints
 PetscErrorCode JacResCopyRes(JacRes *jr, Vec f);
diff --git a/src/LaMEMLib.cpp b/src/LaMEMLib.cpp
index d30e6636..589dd706 100644
--- a/src/LaMEMLib.cpp
+++ b/src/LaMEMLib.cpp
@@ -629,7 +629,7 @@ PetscErrorCode LaMEMLibSolve(LaMEMLib *lm, void *param, PetscLogStage stages[4])
 	PetscCall(PetscLogStagePush(stages[0])); /* Start profiling stage*/
 
 	ierr = LaMEMLibInitGuess(lm, snes); CHKERRQ(ierr);
-    
+
 	PetscCall(PetscLogStagePop()); /* Stop profiling stage*/
 
 	if (param)
@@ -822,10 +822,10 @@ PetscErrorCode LaMEMLibInitGuess(LaMEMLib *lm, SNES snes)
 	ierr = LaMEMLibDiffuseTemp(lm); CHKERRQ(ierr);
 
 	// initialize pressure
-	ierr = JacResInitPres(&lm->jr); CHKERRQ(ierr);
+	ierr = JacResInitPres(&lm->jr,&lm->ts); CHKERRQ(ierr);
 
 	// lithostatic pressure initializtaion
-	ierr = JacResInitLithPres(&lm->jr, &lm->actx); CHKERRQ(ierr);
+	ierr = JacResInitLithPres(&lm->jr, &lm->actx, &lm->ts); CHKERRQ(ierr);
 
 	// compute inverse elastic parameters (dependent on dt)
 	ierr = JacResGetI2Gdt(&lm->jr); CHKERRQ(ierr);
@@ -864,6 +864,7 @@ PetscErrorCode LaMEMLibInitGuess(LaMEMLib *lm, SNES snes)
 PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm)
 {
 	JacRes         *jr;
+	TSSol          *ts;
 	Controls       *ctrl;
 	AdvCtx         *actx;
 	PetscLogDouble t;
@@ -874,12 +875,13 @@ PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm)
 	PetscFunctionBeginUser;
 
 	// access context
+	ts       = &lm->ts; 
 	jr      = &lm->jr;
 	ctrl    = &jr->ctrl;
 	actx    = &lm->actx;
 
 	// check for infinite diffusion
-	if (ctrl->actTemp && ctrl->actSteadyTemp)
+	if (ctrl->actTemp && ctrl->actSteadyTemp && ts->istep==0)
 	{
 		PrintStart(&t,"Computing steady-state temperature distribution", NULL);
 
@@ -903,7 +905,7 @@ PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm)
 	}
 
 	// check for additional limited diffusion
-	if (ctrl->actTemp && ctrl->steadyTempStep)
+	if (ctrl->actTemp && ctrl->steadyTempStep && ts->istep==0)
 	{
 		PrintStart(&t,"Diffusing temperature", NULL);