diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.R b/models/multiphase/d3q27_pf_velocity/Dynamics.R
index b7133e5aa..2023981b4 100644
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.R
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.R
@@ -124,7 +124,7 @@ if (Options$thermo){
 	AddStage("BaseInit" , "Init_distributions", save=Fields$group %in% save_initial)
 	AddStage("calcPhase", "calcPhaseF", save=Fields$name=="PhaseF", load=DensityAll$group %in% load_phase)
 	AddStage("BaseIter" , "Run", save=Fields$group %in% save_iteration, load=DensityAll$group %in% load_iteration )
-	AddStage(name="InitFromFieldsStage", load.densities=TRUE, save.fields=TRUE)
+	AddStage(name="InitFromFieldsStage", load=DensityAll$group %in% "init",read=FALSE, save=Fields$group %in% save_initial_PF)
 	# STAGES FOR VARIOUS OPTIONS
 	if (Options$geometric){
 		AddStage("WallInit_CA"  , "Init_wallNorm", save=Fields$group %in% c("nw", "solid_boundary", extra_fields_to_load_for_bc))
diff --git a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
index 3c18c9fef..02cf4ce43 100755
--- a/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
+++ b/models/multiphase/d3q27_pf_velocity/Dynamics.c.Rt
@@ -284,6 +284,7 @@ CudaDeviceFunction void InitFromFieldsStage()
 	V = Init_UY_External;
 	W = Init_UZ_External;
 	if ( IamWall || IamSolid ) PhaseF = -999;
+	pnorm = 0.0; // initialise as zero and fill in later stage
 }
 
 CudaDeviceFunction void specialCases_Init()