diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index d7aff78d25324f36f466294a3b357b88df9774b4..846cb953b780368598eb7c5a70ed9cfff337711a 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -616,16 +616,18 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	    p_outarray[myompthread] = p_output_2;
 	  } else {
 	    // Thread 0 of non-zero MPI processes needs to allocate memory for the
-	    // output of all threads.
-	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
+	    // output of all threads _doing something_.
+	    int iterstodo = cid_2->number_of_scales - myjxi488 + 1;
+	    if (iterstodo > ompnumthreads) iterstodo = ompnumthreads;
+	    p_output_2 = new InclusionOutputInfo(sconf, gconf, mpidata, myjxi488, iterstodo);
 	    p_outarray[0] = p_output_2;
 	  }
 	  int jer = inclusion_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
 	} else {
-	  if (myompthread > 0) {
+	  //if (myompthread > 0) {
 	    // If there is no input for this thread, set the output pointer to NULL.
 	    p_outarray[myompthread] = NULL;
-	  }	  
+	    //}	  
 	}
 
 #pragma omp barrier
@@ -643,8 +645,18 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  // thread 0 sends the collected virtualfiles to thread 0 of MPI process 0, then deletes them
 	  for (int rr=1; rr<mpidata->nprocs; rr++) {
 	    if (rr == mpidata->rank) {
-	      p_outarray[0]->mpisend(mpidata);
-	      delete p_outarray[0];
+	      if (p_outarray[0] == NULL) {
+		// signal that we are not sending anything
+		int skip_flag = 1;
+		MPI_Send(&skip_flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+	      }
+	      else {
+		// signal that we are sending something
+		int skip_flag = 0;
+		MPI_Send(&skip_flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+		p_outarray[0]->mpisend(mpidata);
+		delete p_outarray[0];
+	      }
 	      vtppoanarray[0]->mpisend(mpidata);
 	      delete vtppoanarray[0];
 	    }
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 3cf5c1b474dcd3eeca1d71677abbd882b54532e9..3f29eb0638bb096bb4a5c24e531e738dba5227ff 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -4485,157 +4485,162 @@ int InclusionOutputInfo::write_legacy(const std::string &output) {
 #ifdef MPI_VERSION
 int InclusionOutputInfo::mpireceive(const mixMPI* mpidata, int pid) {
   int result = 0;
+  int skip_flag;
   int chk_nsph, chk_inpol, chk_iavm, chk_isam, chk_num_theta, chk_num_thetas;
   int chk_num_phi, chk_num_phis, chk_ndirs, chk_idfc, chk_configs;
   double chk_exri;
-  MPI_Recv(&chk_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_iavm, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  result += (chk_nsph == nsph) ? 0 : 1;
-  result += (chk_inpol == inpol) ? 0 : 1;
-  result += (chk_iavm == iavm) ? 0 : 1;
-  result += (chk_isam == isam) ? 0 : 1;
-  result += (chk_num_theta == _num_theta) ? 0 : 1;
-  result += (chk_num_thetas == _num_thetas) ? 0 : 1;
-  result += (chk_num_phi == _num_phi) ? 0 : 1;
-  result += (chk_num_phis == _num_phis) ? 0 : 1;
-  result += (chk_ndirs == ndirs) ? 0 : 1;
-  result += (chk_exri == exri) ? 0 : 1;
-  result += (chk_idfc == idfc) ? 0 : 1;
-  result += (chk_configs == configurations) ? 0 : 1;
-  if (result == 0) {
-    int xi1, offset, chunk_size;
-    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
-    MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    // Receive vectors of single values per scale
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = xi1 - _first_xi;
-    MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_raprs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_raprs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  MPI_Recv(&skip_flag, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  // Proceed with the rest _only if__ skip_flag==0, else nothing is to be received
+  if (skip_flag == 0) {
+    MPI_Recv(&chk_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_iavm, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    result += (chk_nsph == nsph) ? 0 : 1;
+    result += (chk_inpol == inpol) ? 0 : 1;
+    result += (chk_iavm == iavm) ? 0 : 1;
+    result += (chk_isam == isam) ? 0 : 1;
+    result += (chk_num_theta == _num_theta) ? 0 : 1;
+    result += (chk_num_thetas == _num_thetas) ? 0 : 1;
+    result += (chk_num_phi == _num_phi) ? 0 : 1;
+    result += (chk_num_phis == _num_phis) ? 0 : 1;
+    result += (chk_ndirs == ndirs) ? 0 : 1;
+    result += (chk_exri == exri) ? 0 : 1;
+    result += (chk_idfc == idfc) ? 0 : 1;
+    result += (chk_configs == configurations) ? 0 : 1;
+    if (result == 0) {
+      int xi1, offset, chunk_size;
+      MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+      MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      // Receive vectors of single values per scale
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = xi1 - _first_xi;
+      MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_raprs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_raprs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
 
-    // Receive vectors whose sizes depend on configurations
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = (xi1 - _first_xi) * (configurations + 1);
-    MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      // Receive vectors whose sizes depend on configurations
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = (xi1 - _first_xi) * (configurations + 1);
+      MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
 
-    // Receive vectors whose sizes depend on directions and wavelengths
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = (xi1 - _first_xi) * ndirs;
-    MPI_Recv(vec_dir_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_rapr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_rapr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      // Receive vectors whose sizes depend on directions and wavelengths
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = (xi1 - _first_xi) * ndirs;
+      MPI_Recv(vec_dir_scs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_scs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_abs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_abs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_exs1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_exs2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_albeds1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_albeds2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_scsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_scsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_absrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_absrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_exsrt1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_exsrt2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fsas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fsas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fsas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fsas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_qschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_qschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_pschu1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_pschu2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_s0mag1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_s0mag2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_cosav1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_cosav2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_rapr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_rapr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
 
-    MPI_Recv(vec_dir_fl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqel1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqel2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqer1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqer2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqex1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqex2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqey1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqey2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqez1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqez2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_tqsz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_mull + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_mulllr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  }  else {
-    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+      MPI_Recv(vec_dir_fl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqel1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqel2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqer1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqer2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqex1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqex2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqey1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqey2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqez1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqez2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsl1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsl2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsr1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsr2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsx1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsx2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsy1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsy2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsz1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_tqsz2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_mull + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_mulllr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    }  else {
+      MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+    }
   }
   return result;
 }
@@ -5741,90 +5746,95 @@ int SphereOutputInfo::write_legacy(const std::string &file_name) {
 #ifdef MPI_VERSION
 int SphereOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
   int result = 0;
+  int skip_flag;
   int chk_nsph, chk_inpol, chk_isam, chk_num_theta, chk_num_thetas;
   int chk_num_phi, chk_num_phis, chk_ndirs, chk_idfc, chk_configs;
   double chk_exri;
-  MPI_Recv(&chk_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  result += (chk_nsph == nsph) ? 0 : 1;
-  result += (chk_inpol == inpol) ? 0 : 1;
-  result += (chk_isam == isam) ? 0 : 1;
-  result += (chk_num_theta == _num_theta) ? 0 : 1;
-  result += (chk_num_thetas == _num_thetas) ? 0 : 1;
-  result += (chk_num_phi == _num_phi) ? 0 : 1;
-  result += (chk_num_phis == _num_phis) ? 0 : 1;
-  result += (chk_ndirs == ndirs) ? 0 : 1;
-  result += (chk_exri == exri) ? 0 : 1;
-  result += (chk_idfc == idfc) ? 0 : 1;
-  result += (chk_configs == configurations) ? 0 : 1;
-  if (result == 0) {
-    int xi1, offset, chunk_size;
-    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
-    MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    // Receive vectors of single values per scale
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = xi1 - _first_xi;
-    MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    if (nsph != 1) {
-      MPI_Recv(vec_fsat + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE); 
-      MPI_Recv(vec_qschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-      MPI_Recv(vec_pschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-      MPI_Recv(vec_s0magt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    }
+  MPI_Recv(&skip_flag, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  // Proceed with the rest _only if__ skip_flag==0, else nothing is to be received
+  if (skip_flag == 0) {
+    MPI_Recv(&chk_nsph, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_inpol, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_isam, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_theta, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_thetas, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_phi, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_num_phis, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_ndirs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_exri, 1, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_idfc, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    MPI_Recv(&chk_configs, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    result += (chk_nsph == nsph) ? 0 : 1;
+    result += (chk_inpol == inpol) ? 0 : 1;
+    result += (chk_isam == isam) ? 0 : 1;
+    result += (chk_num_theta == _num_theta) ? 0 : 1;
+    result += (chk_num_thetas == _num_thetas) ? 0 : 1;
+    result += (chk_num_phi == _num_phi) ? 0 : 1;
+    result += (chk_num_phis == _num_phis) ? 0 : 1;
+    result += (chk_ndirs == ndirs) ? 0 : 1;
+    result += (chk_exri == exri) ? 0 : 1;
+    result += (chk_idfc == idfc) ? 0 : 1;
+    result += (chk_configs == configurations) ? 0 : 1;
+    if (result == 0) {
+      int xi1, offset, chunk_size;
+      MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+      MPI_Recv(&xi1, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      // Receive vectors of single values per scale
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = xi1 - _first_xi;
+      MPI_Recv(vec_jxi + offset, chunk_size, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_ier + offset, chunk_size, MPI_SHORT, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_vk + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_xi + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      if (nsph != 1) {
+	MPI_Recv(vec_fsat + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE); 
+	MPI_Recv(vec_qschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+	MPI_Recv(vec_pschut + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+	MPI_Recv(vec_s0magt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      }
 
-    // Receive vectors whose sizes depend on configurations and wavelengths
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = (xi1 - _first_xi) * configurations;
-    MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_abs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_albeds + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_scsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_absrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_exsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_fsas + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_qschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_pschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_s0mag + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_cosav + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_raprs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-      
-    // Receive vectors whose sizes depend on NSPH, directions and wavelengths
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = (xi1 - _first_xi) * nsph * ndirs;
-    MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_muls + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_mulslr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      // Receive vectors whose sizes depend on configurations and wavelengths
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = (xi1 - _first_xi) * configurations;
+      MPI_Recv(vec_sphere_sizes + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_sphere_ref_indices + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_abs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_albeds + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_scsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_absrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_exsrt + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_fsas + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_qschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_pschu + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_s0mag + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_cosav + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_raprs + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_tqek1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_tqek2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_tqsk1 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_tqsk2 + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
 
-    // Receive vectors whose sizes depend on NSPH, incidence directions and wavelengths
-    MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
-    MPI_Recv(vec_dir_fx + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fy + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-    MPI_Recv(vec_dir_fz + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  }  else {
-    MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+      // Receive vectors whose sizes depend on NSPH, directions and wavelengths
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = (xi1 - _first_xi) * nsph * ndirs;
+      MPI_Recv(vec_dir_sas11 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas21 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas12 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_sas22 + offset, chunk_size, MPI_C_DOUBLE_COMPLEX, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_muls + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_mulslr + 16 * offset, 16 * chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+      // Receive vectors whose sizes depend on NSPH, incidence directions and wavelengths
+      MPI_Recv(&chunk_size, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
+      MPI_Recv(vec_dir_fx + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fy + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      MPI_Recv(vec_dir_fz + offset, chunk_size, MPI_DOUBLE, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    }  else {
+      MPI_Send(&result, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD);
+    }
   }
   return result;
 }
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index a2d5228ec09e24aa32772daf739bf0785791e836..18d19b80bc716356e40f39449d21bde53c2af666 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -484,16 +484,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou
 	    p_outarray[myompthread] = p_output_2;
 	  } else {
 	    // Thread 0 of non-zero MPI processes needs to allocate memory for the
-	    // output of all threads.
-	    p_output_2 = new SphereOutputInfo(sconf, gconf, mpidata, myjxi488, ompnumthreads);
+	    // output of all threads _doing something_.
+	    int iterstodo = sid_2->number_of_scales - myjxi488 + 1;
+	    if (iterstodo > ompnumthreads) iterstodo = ompnumthreads;
+	    p_output_2 = new SphereOutputInfo(sconf, gconf, mpidata, myjxi488, iterstodo);
 	    p_outarray[0] = p_output_2;
 	  }
 	  int jer = sphere_jxi488_cycle(myjxi488 - 1, sconf, gconf, p_sa, sid_2, p_output_2, output_path, vtppoanp_2);
 	} else {
-	  if (myompthread > 0) {
+	  // if (myompthread > 0) {
 	    // If there is no input for this thread, set the output pointer to NULL.
 	    p_outarray[myompthread] = NULL;
-	  }	  
+	    //}	  
 	}
 
 #pragma omp barrier
@@ -511,8 +513,18 @@ void sphere(const string& config_file, const string& data_file, const string& ou
 	  // thread 0 sends the collected virtualfiles to thread 0 of MPI process 0, then deletes them
 	  for (int rr=1; rr<mpidata->nprocs; rr++) {
 	    if (rr == mpidata->rank) {
-	      p_outarray[0]->mpisend(mpidata);
-	      delete p_outarray[0];
+	      if (p_outarray[0] == NULL) {
+		// signal that we are not sending anything
+		int skip_flag = 1;
+		MPI_Send(&skip_flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+	      }
+	      else {
+		// signal that we are sending something
+		int skip_flag = 0;
+		MPI_Send(&skip_flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+		p_outarray[0]->mpisend(mpidata);
+		delete p_outarray[0];
+	      }
 	      vtppoanarray[0]->mpisend(mpidata);
 	      delete vtppoanarray[0];
 	    }