diff --git a/spm_mb_fit.m b/spm_mb_fit.m
index 38847a1..667ab22 100644
--- a/spm_mb_fit.m
+++ b/spm_mb_fit.m
@@ -17,15 +17,6 @@
 %--------------------------------------------------------------------------
 rng('default'); rng(1);
 
-% If SPM has been compiled with OpenMP support then the number of threads
-% are here set to speed up the algorithm
-%--------------------------------------------------------------------------
-if sett.nworker > 1
-    setenv('SPM_NUM_THREADS',sprintf('%d',0));
-else
-    setenv('SPM_NUM_THREADS',sprintf('%d',-1));
-end
-
 % Get template size and orientation
 %--------------------------------------------------------------------------
 if isfield(sett.mu,'exist')
@@ -40,6 +31,24 @@
 sz      = spm_mb_shape('zoom_settings',sett.v_settings,sett.mu,nz);
 sett.ms = sz(end);
 
+% Configure number of parfor workers
+%--------------------------------------------------------------------------
+nworker_zm = sett.nworker;
+if isscalar(nworker_zm), nworker_zm = [nworker_zm]; end
+if numel(nworker_zm) > numel(sz), nworker_zm = nworker_zm(1:numel(sz)); end
+nworker_zm = padarray(nworker_zm,[0, numel(sz) - numel(nworker_zm)],nworker_zm(1),'pre');
+nworker_zm = fliplr(nworker_zm);   
+sett.nworker = nworker_zm(end);    
+    
+% If SPM has been compiled with OpenMP support then the number of threads
+% are here set to speed up the algorithm
+%--------------------------------------------------------------------------
+if sett.nworker > 1
+    setenv('SPM_NUM_THREADS',sprintf('%d',0));
+else
+    setenv('SPM_NUM_THREADS',sprintf('%d',-1));
+end
+
 % Init shape model parameters
 %--------------------------------------------------------------------------
 dat = spm_mb_shape('init_def',dat,sett.ms);
@@ -78,6 +87,7 @@
 spm_plot_convergence('Init','Rigid Alignment','Objective','Iteration');
 E      = Inf;
 for it0=1:nit_aff
+   
     if it0>1
         oE  = E/nvox(dat);
     else
@@ -125,6 +135,9 @@
 for zm=numel(sz):-1:1 % loop over zoom levels
     fprintf('\nzoom=1/%d: %d x %d x %d\n', 2^(zm-1), sett.ms.d);
 
+    sett.nworker = nworker_zm(zm);
+    disp(['sett.nworker = ' num2str(sett.nworker)])
+    
     if updt_mu
         dat = spm_mb_appearance('restart',dat,sett);
     end
@@ -202,6 +215,8 @@
     end
     do_save(mu,sett,dat);
 end
+% reset number of workers
+sett.nworker = nworker_zm(end);
 %spm_plot_convergence('Clear');
 %==========================================================================
 
diff --git a/spm_mb_output.m b/spm_mb_output.m
index f2ebef8..d34b5e7 100644
--- a/spm_mb_output.m
+++ b/spm_mb_output.m
@@ -24,6 +24,8 @@
 % If SPM has been compiled with OpenMP support then the number of threads
 % are here set to speed up the algorithm
 %--------------------------------------------------------------------------
+if isscalar(sett.nworker), sett.nworker = [sett.nworker]; end
+sett.nworker = sett.nworker(end);
 nw   = spm_mb_shape('get_num_workers',sett,max(27,sett.K*5+17));
 if sett.nworker > 1
     setenv('SPM_NUM_THREADS',sprintf('%d',0));