diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 46de0e4970fd32402cd1b9b96d1b9ae902131263..f3e9141dd5567e96a0155c03099870a1f6d7fcd8 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -440,8 +440,8 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	    int jer = cluster_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
 	  } else {
 	    if (myompthread > 0) {
-	      // If there is no input for this thread, set output pointer to NULL.
-	      p_outarray[myompthread] = NULL;
+	      // If there is no input for this thread, mark to skip.
+	      p_outarray[myompthread] = new ClusterOutputInfo(1);
 	    }
 	  }
 #pragma omp barrier
@@ -453,11 +453,9 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	  // threads different from 0 append their virtual files to the one of thread 0, and delete them
 	  if (myompthread == 0) {
 	    for (int ti=1; ti<ompnumthreads; ti++) {
-	      if (p_outarray[ti] != NULL) {
-		p_outarray[0]->insert(*(p_outarray[ti]));
-		delete p_outarray[ti];
-		p_outarray[ti] = NULL;
-	      }
+	      p_outarray[0]->insert(*(p_outarray[ti]));
+	      delete p_outarray[ti];
+	      p_outarray[ti] = NULL;
 	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	      delete vtppoanarray[ti];
 	    }
@@ -628,39 +626,25 @@ void cluster(const string& config_file, const string& data_file, const string& o
 	  }
 	  int jer = cluster_jxi488_cycle(myjxi488, sconf, gconf, p_scattering_angles, cid_2, p_output_2, output_path, vtppoanp_2);
 	} else {
-	  // if (myompthread > 0) {
-	    // If there is no input for this thread, set the output pointer to NULL.
-	    p_outarray[myompthread] = NULL;
-	  //}
+	  p_outarray[myompthread] = new ClusterOutputInfo(1);
 	}
 
 #pragma omp barrier
 	// threads different from 0 append their virtual files to the one of thread 0, and delete them
 	if (myompthread == 0) {
 	  for (int ti=1; ti<ompnumthreads; ti++) {
-	    if (p_outarray[ti] != NULL) {
-	      p_outarray[0]->insert(*(p_outarray[ti]));
-	      delete p_outarray[ti];
-	      p_outarray[ti] = NULL;
-	    }
+	    p_outarray[0]->insert(*(p_outarray[ti]));
+	    delete p_outarray[ti];
+	    p_outarray[ti] = NULL;
 	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	    delete vtppoanarray[ti];
 	  }
 	  // 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) {
-	      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];
-	      }
+	      p_outarray[0]->mpisend(mpidata);
+	      delete p_outarray[0];
+	      p_outarray[0] = NULL;
 	      vtppoanarray[0]->mpisend(mpidata);
 	      delete vtppoanarray[0];
 	    }
diff --git a/src/include/errors.h b/src/include/errors.h
index a986defa166fb9c7edbcdf5911f45084461e6590..1908ade318ecb90db059bde369971a27bcf85d16 100644
--- a/src/include/errors.h
+++ b/src/include/errors.h
@@ -22,6 +22,32 @@
 #ifndef INCLUDE_ERRORS_H_
 #define INCLUDE_ERRORS_H_
 
+/*! \brief Exception for wrong OutputInfo NULL calls.
+ */
+class UnrecognizedOutputInfo: public std::exception {
+protected:
+  //! Description of the problem.
+  std::string message;
+
+public:
+  /**
+   * \brief Exception instance constructor.
+   *
+   * \param requested: `int` The index that was requested.
+   */
+  UnrecognizedOutputInfo(int requested) {
+    message = "Error: passed parameter " + std::to_string(requested)
+      + ", but only 1 is allowed";
+  }
+  
+  /**
+   * \brief Exception message.
+   */
+  virtual const char* what() const throw() {
+    return message.c_str();
+  }
+};
+
 /*! \brief Exception for out of bounds List requests.
  */
 class ListOutOfBoundsException: public std::exception {
diff --git a/src/include/outputs.h b/src/include/outputs.h
index b4f9b02ed34cbb91a56cfbe9de2832e6fd28be52..82ee40ae0c2d19dbaa18135e680c1f001a94ecd6 100644
--- a/src/include/outputs.h
+++ b/src/include/outputs.h
@@ -33,6 +33,8 @@
  */
 class ClusterOutputInfo {
 protected:
+  //! \brief Flag for skipping mpisend() and mpireceive()
+  int _skip_flag;
   //! \brief Number of incident azimuth calculations.
   int _num_theta;
   //! \brief Number of scattered azimuth calculations.
@@ -63,6 +65,8 @@ protected:
   int write_legacy(const std::string &output);
   
 public:
+  //! \brief Read-only view on skip_flag
+  const int &skip_flag = _skip_flag;
   //! \brief Read-only view on the ID of the first scale
   const int &first_xi = _first_xi;
   //! \brief Number of spheres in the aggregate.
@@ -471,6 +475,12 @@ public:
    */   
   ClusterOutputInfo(const std::string &hdf5_name);
 
+  /*! \brief `ClusterOutputInfo` constructor for the dummy NULL case
+   *
+   * \param skip_flag: `const int` must be passed as the `1` constant.
+   */   
+  ClusterOutputInfo(const int skip_flag);
+
   /*! \brief `ClusterOutputInfo` instance destroyer.
    */
   ~ClusterOutputInfo();
@@ -551,6 +561,8 @@ public:
  */
 class InclusionOutputInfo {
 protected:
+  //! \brief Flag for skipping mpisend() and mpireceive()
+  int _skip_flag;
   //! \brief Number of incident azimuth calculations.
   int _num_theta;
   //! \brief Number of scattered azimuth calculations.
@@ -581,6 +593,8 @@ protected:
   int write_legacy(const std::string &output);
   
 public:
+  //! \brief Read-only view on skip_flag
+  const int &skip_flag = _skip_flag;
   //! \brief Read-only view on the ID of the first scale
   const int &first_xi = _first_xi;
   //! \brief Number of spheres in the aggregate.
@@ -901,6 +915,12 @@ public:
    */   
   InclusionOutputInfo(const std::string &hdf5_name);
 
+  /*! \brief `InclusionOutputInfo` constructor for the dummy NULL case
+   *
+   * \param skip_flag: `const int` must be passed as the `1` constant.
+   */   
+  InclusionOutputInfo(const int skip_flag);
+
   /*! \brief `InclusionOutputInfo` instance destroyer.
    */
   ~InclusionOutputInfo();
@@ -981,6 +1001,8 @@ public:
  */
 class SphereOutputInfo {
 protected:
+  //! \brief Flag for skipping mpisend() and mpireceive()
+  int _skip_flag;
   //! \brief Number of incident azimuth calculations.
   int _num_theta;
   //! \brief Number of scattered azimuth calculations.
@@ -1011,6 +1033,8 @@ protected:
   int write_legacy(const std::string &output);
   
 public:
+  //! \brief Read-only view on skip_flag
+  const int &skip_flag = _skip_flag;
   //! \brief Read-only view on the ID of the first scale
   const int &first_xi = _first_xi;
   //! \brief Number of spheres.
@@ -1181,6 +1205,12 @@ public:
    */   
   SphereOutputInfo(const std::string &hdf5_name);
 
+  /*! \brief `SphereOutputInfo` constructor for the dummy NULL case
+   *
+   * \param skip_flag: `const int` must be passed as the `1` constant.
+   */   
+  SphereOutputInfo(const int skip_flag);
+
   /*! \brief `InclusionOutputInfo` instance destroyer.
    */
   ~SphereOutputInfo();
diff --git a/src/inclusion/inclusion.cpp b/src/inclusion/inclusion.cpp
index 846cb953b780368598eb7c5a70ed9cfff337711a..0b0033c4a20f7609c5fd5e9dc940efdf8f7509df 100644
--- a/src/inclusion/inclusion.cpp
+++ b/src/inclusion/inclusion.cpp
@@ -437,7 +437,7 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  } else {
 	    if (myompthread > 0) {
 	      // If there is no input for this thread, set output pointer to NULL.
-	      p_outarray[myompthread] = NULL;
+	      p_outarray[myompthread] = new InclusionOutputInfo(1);
 	    }
 	  }
 #pragma omp barrier
@@ -449,11 +449,9 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  // threads different from 0 append their virtual files to the one of thread 0, and delete them
 	  if (myompthread == 0) {
 	    for (int ti=1; ti<ompnumthreads; ti++) {
-	      if (p_outarray[ti] != NULL) {
-		p_outarray[0]->insert(*(p_outarray[ti]));
-		delete p_outarray[ti];
-		p_outarray[ti] = NULL;
-	      }
+	      p_outarray[0]->insert(*(p_outarray[ti]));
+	      delete p_outarray[ti];
+	      p_outarray[ti] = NULL;
 	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	      delete vtppoanarray[ti];
 	    }
@@ -624,39 +622,25 @@ void inclusion(const string& config_file, const string& data_file, const string&
 	  }
 	  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 there is no input for this thread, set the output pointer to NULL.
-	    p_outarray[myompthread] = NULL;
-	    //}	  
+	  p_outarray[myompthread] = new InclusionOutputInfo(1);
 	}
 
 #pragma omp barrier
 	// threads different from 0 append their virtual files to the one of thread 0, and delete them
 	if (myompthread == 0) {
 	  for (int ti=1; ti<ompnumthreads; ti++) {
-	    if (p_outarray[ti] != NULL) {
-	      p_outarray[0]->insert(*(p_outarray[ti]));
-	      delete p_outarray[ti];
-	      p_outarray[ti] = NULL;
-	    }
+	    p_outarray[0]->insert(*(p_outarray[ti]));
+	    delete p_outarray[ti];
+	    p_outarray[ti] = NULL;
 	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	    delete vtppoanarray[ti];
 	  }
 	  // 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) {
-	      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];
-	      }
+	      p_outarray[0]->mpisend(mpidata);
+	      delete p_outarray[0];
+	      p_outarray[0] = NULL;
 	      vtppoanarray[0]->mpisend(mpidata);
 	      delete vtppoanarray[0];
 	    }
diff --git a/src/libnptm/outputs.cpp b/src/libnptm/outputs.cpp
index 3f29eb0638bb096bb4a5c24e531e738dba5227ff..48791dc9a947ef24f9183d974b5855d7ebe24602 100644
--- a/src/libnptm/outputs.cpp
+++ b/src/libnptm/outputs.cpp
@@ -61,6 +61,7 @@ ClusterOutputInfo::ClusterOutputInfo(
   ScattererConfiguration *sc, GeometryConfiguration *gc,
   const mixMPI *mpidata, int first_xi, int xi_length
 ) {
+  _skip_flag = 0;
   nsph = gc->number_of_spheres;
   li = gc->li;
   le = gc->le;
@@ -293,6 +294,7 @@ ClusterOutputInfo::ClusterOutputInfo(const std::string &hdf5_name) {
   HDFFile *hdf_file = new HDFFile(hdf5_name, flags);
   herr_t status = hdf_file->get_status();
   string str_name, str_type;
+  _skip_flag = 0;
   if (status == 0) {
     status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
     status = hdf_file->read("LI", "INT32_(1)", &li);
@@ -701,170 +703,184 @@ ClusterOutputInfo::ClusterOutputInfo(const std::string &hdf5_name) {
   }
 }
 
+ClusterOutputInfo::ClusterOutputInfo(const int flag) {
+  /*
+    create a dummy placeholder just to know I should skip MPI_Send and MPI_Recv
+  */
+  if (flag == 1) {
+    _skip_flag = 1;
+  } else {
+    UnrecognizedOutputInfo ex(flag);
+    throw ex;
+  }
+}
+
 ClusterOutputInfo::~ClusterOutputInfo() {
-  delete[] vec_x_coords;
-  delete[] vec_y_coords;
-  delete[] vec_z_coords;
-  delete[] vec_jxi;
-  delete[] vec_ier;
-  delete[] vec_vk;
-  delete[] vec_xi;
-  delete[] vec_sphere_sizes;
-  delete[] vec_sphere_ref_indices;
-  delete[] vec_sphere_scs;
-  delete[] vec_sphere_abs;
-  delete[] vec_sphere_exs;
-  delete[] vec_sphere_albs;
-  delete[] vec_sphere_sqscs;
-  delete[] vec_sphere_sqabs;
-  delete[] vec_sphere_sqexs;
-  delete[] vec_fsas;
-  delete[] vec_qschus;
-  delete[] vec_pschus;
-  delete[] vec_s0mags;
-  delete[] vec_cosavs;
-  delete[] vec_raprs;
-  delete[] vec_tqek1;
-  delete[] vec_tqsk1;
-  delete[] vec_tqek2;
-  delete[] vec_tqsk2;
-  delete[] vec_fsat;
-  delete[] vec_qschut;
-  delete[] vec_pschut;
-  delete[] vec_s0magt;
-  delete[] vec_scc1;
-  delete[] vec_scc2;
-  delete[] vec_abc1;
-  delete[] vec_abc2;
-  delete[] vec_exc1;
-  delete[] vec_exc2;
-  delete[] vec_albedc1;
-  delete[] vec_albedc2;
-  delete[] vec_sccrt1;
-  delete[] vec_sccrt2;
-  delete[] vec_abcrt1;
-  delete[] vec_abcrt2;
-  delete[] vec_excrt1;
-  delete[] vec_excrt2;
-  delete[] vec_fsac11;
-  delete[] vec_fsac21;
-  delete[] vec_fsac22;
-  delete[] vec_fsac12;
-  delete[] vec_qschuc1;
-  delete[] vec_qschuc2;
-  delete[] vec_pschuc1;
-  delete[] vec_pschuc2;
-  delete[] vec_s0magc1;
-  delete[] vec_s0magc2;
-  delete[] vec_cosavc1;
-  delete[] vec_cosavc2;
-  delete[] vec_raprc1;
-  delete[] vec_raprc2;
-  delete[] vec_fkc1;
-  delete[] vec_fkc2;
-  delete[] vec_dir_tidg;
-  delete[] vec_dir_pidg;
-  delete[] vec_dir_tsdg;
-  delete[] vec_dir_psdg;
-  delete[] vec_dir_scand;
-  delete[] vec_dir_cfmp;
-  delete[] vec_dir_sfmp;
-  delete[] vec_dir_cfsp;
-  delete[] vec_dir_sfsp;
-  delete[] vec_dir_un;
-  delete[] vec_dir_uns;
-  delete[] vec_dir_sas11;
-  delete[] vec_dir_sas21;
-  delete[] vec_dir_sas12;
-  delete[] vec_dir_sas22;
-  delete[] vec_dir_muls;
-  delete[] vec_dir_mulslr;
-  delete[] vec_dir_sat11;
-  delete[] vec_dir_sat21;
-  delete[] vec_dir_sat12;
-  delete[] vec_dir_sat22;
-  delete[] vec_dir_scc1;
-  delete[] vec_dir_scc2;
-  delete[] vec_dir_abc1;
-  delete[] vec_dir_abc2;
-  delete[] vec_dir_exc1;
-  delete[] vec_dir_exc2;
-  delete[] vec_dir_albedc1;
-  delete[] vec_dir_albedc2;
-  delete[] vec_dir_qscc1;
-  delete[] vec_dir_qscc2;
-  delete[] vec_dir_qabc1;
-  delete[] vec_dir_qabc2;
-  delete[] vec_dir_qexc1;
-  delete[] vec_dir_qexc2;
-  delete[] vec_qscamc1;
-  delete[] vec_qscamc2;
-  delete[] vec_qabsmc1;
-  delete[] vec_qabsmc2;
-  delete[] vec_qextmc1;
-  delete[] vec_qextmc2;
-  delete[] vec_dir_sccrt1;
-  delete[] vec_dir_sccrt2;
-  delete[] vec_dir_abcrt1;
-  delete[] vec_dir_abcrt2;
-  delete[] vec_dir_excrt1;
-  delete[] vec_dir_excrt2;
-  delete[] vec_dir_fsac11;
-  delete[] vec_dir_fsac21;
-  delete[] vec_dir_fsac12;
-  delete[] vec_dir_fsac22;
-  delete[] vec_dir_sac11;
-  delete[] vec_dir_sac21;
-  delete[] vec_dir_sac12;
-  delete[] vec_dir_sac22;
-  delete[] vec_dir_qschuc1;
-  delete[] vec_dir_qschuc2;
-  delete[] vec_dir_pschuc1;
-  delete[] vec_dir_pschuc2;
-  delete[] vec_dir_s0magc1;
-  delete[] vec_dir_s0magc2;
-  delete[] vec_dir_cosavc1;
-  delete[] vec_dir_cosavc2;
-  delete[] vec_dir_raprc1;
-  delete[] vec_dir_raprc2;
-  delete[] vec_dir_flc1;
-  delete[] vec_dir_flc2;
-  delete[] vec_dir_frc1;
-  delete[] vec_dir_frc2;
-  delete[] vec_dir_fkc1;
-  delete[] vec_dir_fkc2;
-  delete[] vec_dir_fxc1;
-  delete[] vec_dir_fxc2;
-  delete[] vec_dir_fyc1;
-  delete[] vec_dir_fyc2;
-  delete[] vec_dir_fzc1;
-  delete[] vec_dir_fzc2;
-  delete[] vec_dir_tqelc1;
-  delete[] vec_dir_tqelc2;
-  delete[] vec_dir_tqerc1;
-  delete[] vec_dir_tqerc2;
-  delete[] vec_dir_tqekc1;
-  delete[] vec_dir_tqekc2;
-  delete[] vec_dir_tqexc1;
-  delete[] vec_dir_tqexc2;
-  delete[] vec_dir_tqeyc1;
-  delete[] vec_dir_tqeyc2;
-  delete[] vec_dir_tqezc1;
-  delete[] vec_dir_tqezc2;
-  delete[] vec_dir_tqslc1;
-  delete[] vec_dir_tqslc2;
-  delete[] vec_dir_tqsrc1;
-  delete[] vec_dir_tqsrc2;
-  delete[] vec_dir_tqskc1;
-  delete[] vec_dir_tqskc2;
-  delete[] vec_dir_tqsxc1;
-  delete[] vec_dir_tqsxc2;
-  delete[] vec_dir_tqsyc1;
-  delete[] vec_dir_tqsyc2;
-  delete[] vec_dir_tqszc1;
-  delete[] vec_dir_tqszc2;
-  delete[] vec_dir_mulc;
-  delete[] vec_dir_mulclr;
+  if (_skip_flag != 1) {
+    delete[] vec_x_coords;
+    delete[] vec_y_coords;
+    delete[] vec_z_coords;
+    delete[] vec_jxi;
+    delete[] vec_ier;
+    delete[] vec_vk;
+    delete[] vec_xi;
+    delete[] vec_sphere_sizes;
+    delete[] vec_sphere_ref_indices;
+    delete[] vec_sphere_scs;
+    delete[] vec_sphere_abs;
+    delete[] vec_sphere_exs;
+    delete[] vec_sphere_albs;
+    delete[] vec_sphere_sqscs;
+    delete[] vec_sphere_sqabs;
+    delete[] vec_sphere_sqexs;
+    delete[] vec_fsas;
+    delete[] vec_qschus;
+    delete[] vec_pschus;
+    delete[] vec_s0mags;
+    delete[] vec_cosavs;
+    delete[] vec_raprs;
+    delete[] vec_tqek1;
+    delete[] vec_tqsk1;
+    delete[] vec_tqek2;
+    delete[] vec_tqsk2;
+    delete[] vec_fsat;
+    delete[] vec_qschut;
+    delete[] vec_pschut;
+    delete[] vec_s0magt;
+    delete[] vec_scc1;
+    delete[] vec_scc2;
+    delete[] vec_abc1;
+    delete[] vec_abc2;
+    delete[] vec_exc1;
+    delete[] vec_exc2;
+    delete[] vec_albedc1;
+    delete[] vec_albedc2;
+    delete[] vec_sccrt1;
+    delete[] vec_sccrt2;
+    delete[] vec_abcrt1;
+    delete[] vec_abcrt2;
+    delete[] vec_excrt1;
+    delete[] vec_excrt2;
+    delete[] vec_fsac11;
+    delete[] vec_fsac21;
+    delete[] vec_fsac22;
+    delete[] vec_fsac12;
+    delete[] vec_qschuc1;
+    delete[] vec_qschuc2;
+    delete[] vec_pschuc1;
+    delete[] vec_pschuc2;
+    delete[] vec_s0magc1;
+    delete[] vec_s0magc2;
+    delete[] vec_cosavc1;
+    delete[] vec_cosavc2;
+    delete[] vec_raprc1;
+    delete[] vec_raprc2;
+    delete[] vec_fkc1;
+    delete[] vec_fkc2;
+    delete[] vec_dir_tidg;
+    delete[] vec_dir_pidg;
+    delete[] vec_dir_tsdg;
+    delete[] vec_dir_psdg;
+    delete[] vec_dir_scand;
+    delete[] vec_dir_cfmp;
+    delete[] vec_dir_sfmp;
+    delete[] vec_dir_cfsp;
+    delete[] vec_dir_sfsp;
+    delete[] vec_dir_un;
+    delete[] vec_dir_uns;
+    delete[] vec_dir_sas11;
+    delete[] vec_dir_sas21;
+    delete[] vec_dir_sas12;
+    delete[] vec_dir_sas22;
+    delete[] vec_dir_muls;
+    delete[] vec_dir_mulslr;
+    delete[] vec_dir_sat11;
+    delete[] vec_dir_sat21;
+    delete[] vec_dir_sat12;
+    delete[] vec_dir_sat22;
+    delete[] vec_dir_scc1;
+    delete[] vec_dir_scc2;
+    delete[] vec_dir_abc1;
+    delete[] vec_dir_abc2;
+    delete[] vec_dir_exc1;
+    delete[] vec_dir_exc2;
+    delete[] vec_dir_albedc1;
+    delete[] vec_dir_albedc2;
+    delete[] vec_dir_qscc1;
+    delete[] vec_dir_qscc2;
+    delete[] vec_dir_qabc1;
+    delete[] vec_dir_qabc2;
+    delete[] vec_dir_qexc1;
+    delete[] vec_dir_qexc2;
+    delete[] vec_qscamc1;
+    delete[] vec_qscamc2;
+    delete[] vec_qabsmc1;
+    delete[] vec_qabsmc2;
+    delete[] vec_qextmc1;
+    delete[] vec_qextmc2;
+    delete[] vec_dir_sccrt1;
+    delete[] vec_dir_sccrt2;
+    delete[] vec_dir_abcrt1;
+    delete[] vec_dir_abcrt2;
+    delete[] vec_dir_excrt1;
+    delete[] vec_dir_excrt2;
+    delete[] vec_dir_fsac11;
+    delete[] vec_dir_fsac21;
+    delete[] vec_dir_fsac12;
+    delete[] vec_dir_fsac22;
+    delete[] vec_dir_sac11;
+    delete[] vec_dir_sac21;
+    delete[] vec_dir_sac12;
+    delete[] vec_dir_sac22;
+    delete[] vec_dir_qschuc1;
+    delete[] vec_dir_qschuc2;
+    delete[] vec_dir_pschuc1;
+    delete[] vec_dir_pschuc2;
+    delete[] vec_dir_s0magc1;
+    delete[] vec_dir_s0magc2;
+    delete[] vec_dir_cosavc1;
+    delete[] vec_dir_cosavc2;
+    delete[] vec_dir_raprc1;
+    delete[] vec_dir_raprc2;
+    delete[] vec_dir_flc1;
+    delete[] vec_dir_flc2;
+    delete[] vec_dir_frc1;
+    delete[] vec_dir_frc2;
+    delete[] vec_dir_fkc1;
+    delete[] vec_dir_fkc2;
+    delete[] vec_dir_fxc1;
+    delete[] vec_dir_fxc2;
+    delete[] vec_dir_fyc1;
+    delete[] vec_dir_fyc2;
+    delete[] vec_dir_fzc1;
+    delete[] vec_dir_fzc2;
+    delete[] vec_dir_tqelc1;
+    delete[] vec_dir_tqelc2;
+    delete[] vec_dir_tqerc1;
+    delete[] vec_dir_tqerc2;
+    delete[] vec_dir_tqekc1;
+    delete[] vec_dir_tqekc2;
+    delete[] vec_dir_tqexc1;
+    delete[] vec_dir_tqexc2;
+    delete[] vec_dir_tqeyc1;
+    delete[] vec_dir_tqeyc2;
+    delete[] vec_dir_tqezc1;
+    delete[] vec_dir_tqezc2;
+    delete[] vec_dir_tqslc1;
+    delete[] vec_dir_tqslc2;
+    delete[] vec_dir_tqsrc1;
+    delete[] vec_dir_tqsrc2;
+    delete[] vec_dir_tqskc1;
+    delete[] vec_dir_tqskc2;
+    delete[] vec_dir_tqsxc1;
+    delete[] vec_dir_tqsxc2;
+    delete[] vec_dir_tqsyc1;
+    delete[] vec_dir_tqsyc2;
+    delete[] vec_dir_tqszc1;
+    delete[] vec_dir_tqszc2;
+    delete[] vec_dir_mulc;
+    delete[] vec_dir_mulclr;
+  }
 }
 
 long ClusterOutputInfo::compute_size(
@@ -961,182 +977,184 @@ long ClusterOutputInfo::compute_size() {
 
 int ClusterOutputInfo::insert(const ClusterOutputInfo &rhs) {
   int result = 0;
-  result += (rhs.nsph == nsph) ? 0 : 1;
-  result += (rhs.inpol == inpol) ? 0 : 1;
-  result += (rhs.iavm == iavm) ? 0 : 1;
-  result += (rhs.isam == isam) ? 0 : 1;
-  result += (rhs._num_theta == _num_theta) ? 0 : 1;
-  result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
-  result += (rhs._num_phi == _num_phi) ? 0 : 1;
-  result += (rhs._num_phis == _num_phis) ? 0 : 1;
-  result += (rhs.ndirs == ndirs) ? 0 : 1;
-  result += (rhs.exri == exri) ? 0 : 1;
-  result += (rhs.idfc == idfc) ? 0 : 1;
-  result += (rhs.configurations == configurations) ? 0 : 1;
-  if (result == 0) {
-    int offset, chunk_size, xi1;
-    xi1 = rhs._first_xi;
-    // Insert vectors whose size depends on wavelengths
-    offset = xi1 - _first_xi;
-    chunk_size = rhs.xi_block_size;
-    memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
-    memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
-    memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
-    memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
-    memcpy(vec_fsat + offset, rhs.vec_fsat, chunk_size * sizeof(dcomplex));
-    memcpy(vec_qschut + offset, rhs.vec_qschut, chunk_size * sizeof(double));
-    memcpy(vec_pschut + offset, rhs.vec_pschut, chunk_size * sizeof(double));
-    memcpy(vec_s0magt + offset, rhs.vec_s0magt, chunk_size * sizeof(double));
-    memcpy(vec_scc1 + offset, rhs.vec_scc1, chunk_size * sizeof(double));
-    memcpy(vec_scc2 + offset, rhs.vec_scc2, chunk_size * sizeof(double));
-    memcpy(vec_abc1 + offset, rhs.vec_abc1, chunk_size * sizeof(double));
-    memcpy(vec_abc2 + offset, rhs.vec_abc2, chunk_size * sizeof(double));
-    memcpy(vec_exc1 + offset, rhs.vec_exc1, chunk_size * sizeof(double));
-    memcpy(vec_exc2 + offset, rhs.vec_exc2, chunk_size * sizeof(double));
-    memcpy(vec_albedc1 + offset, rhs.vec_albedc1, chunk_size * sizeof(double));
-    memcpy(vec_albedc2 + offset, rhs.vec_albedc2, chunk_size * sizeof(double));
-    memcpy(vec_qscamc1 + offset, rhs.vec_qscamc1, chunk_size * sizeof(double));
-    memcpy(vec_qscamc2 + offset, rhs.vec_qscamc2, chunk_size * sizeof(double));
-    memcpy(vec_qabsmc1 + offset, rhs.vec_qabsmc1, chunk_size * sizeof(double));
-    memcpy(vec_qabsmc2 + offset, rhs.vec_qabsmc2, chunk_size * sizeof(double));
-    memcpy(vec_qextmc1 + offset, rhs.vec_qextmc1, chunk_size * sizeof(double));
-    memcpy(vec_qextmc2 + offset, rhs.vec_qextmc2, chunk_size * sizeof(double));
-    memcpy(vec_sccrt1 + offset, rhs.vec_sccrt1, chunk_size * sizeof(double));
-    memcpy(vec_sccrt2 + offset, rhs.vec_sccrt2, chunk_size * sizeof(double));
-    memcpy(vec_abcrt1 + offset, rhs.vec_abcrt1, chunk_size * sizeof(double));
-    memcpy(vec_abcrt2 + offset, rhs.vec_abcrt2, chunk_size * sizeof(double));
-    memcpy(vec_excrt1 + offset, rhs.vec_excrt1, chunk_size * sizeof(double));
-    memcpy(vec_excrt2 + offset, rhs.vec_excrt2, chunk_size * sizeof(double));
-    memcpy(vec_fsac11 + offset, rhs.vec_fsac11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsac21 + offset, rhs.vec_fsac21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsac22 + offset, rhs.vec_fsac22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsac12 + offset, rhs.vec_fsac12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_qschuc1 + offset, rhs.vec_qschuc1, chunk_size * sizeof(double));
-    memcpy(vec_qschuc2 + offset, rhs.vec_qschuc2, chunk_size * sizeof(double));
-    memcpy(vec_pschuc1 + offset, rhs.vec_pschuc1, chunk_size * sizeof(double));
-    memcpy(vec_pschuc2 + offset, rhs.vec_pschuc2, chunk_size * sizeof(double));
-    memcpy(vec_s0magc1 + offset, rhs.vec_s0magc1, chunk_size * sizeof(double));
-    memcpy(vec_s0magc2 + offset, rhs.vec_s0magc2, chunk_size * sizeof(double));
-    memcpy(vec_cosavc1 + offset, rhs.vec_cosavc1, chunk_size * sizeof(double));
-    memcpy(vec_cosavc2 + offset, rhs.vec_cosavc2, chunk_size * sizeof(double));
-    memcpy(vec_raprc1 + offset, rhs.vec_raprc1, chunk_size * sizeof(double));
-    memcpy(vec_raprc2 + offset, rhs.vec_raprc2, chunk_size * sizeof(double));
-    memcpy(vec_fkc1 + offset, rhs.vec_fkc1, chunk_size * sizeof(double));
-    memcpy(vec_fkc2 + offset, rhs.vec_fkc2, chunk_size * sizeof(double));
-    // Insert vectors of multiple configuration values per scale
-    offset = (xi1 - _first_xi) * configurations;
-    chunk_size = rhs.xi_block_size * configurations;
-    memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
-    memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
-    memcpy(vec_sphere_scs + offset, rhs.vec_sphere_scs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_abs + offset, rhs.vec_sphere_abs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_exs + offset, rhs.vec_sphere_exs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_albs + offset, rhs.vec_sphere_albs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_sqscs + offset, rhs.vec_sphere_sqscs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_sqabs + offset, rhs.vec_sphere_sqabs, chunk_size * sizeof(double));
-    memcpy(vec_sphere_sqexs + offset, rhs.vec_sphere_sqexs, chunk_size * sizeof(double));
-    memcpy(vec_fsas + offset, rhs.vec_fsas, chunk_size * sizeof(dcomplex));
-    memcpy(vec_qschus + offset, rhs.vec_qschus, chunk_size * sizeof(double));
-    memcpy(vec_pschus + offset, rhs.vec_pschus, chunk_size * sizeof(double));
-    memcpy(vec_s0mags + offset, rhs.vec_s0mags, chunk_size * sizeof(double));
-    memcpy(vec_cosavs + offset, rhs.vec_cosavs, chunk_size * sizeof(double));
-    memcpy(vec_raprs + offset, rhs.vec_raprs, chunk_size * sizeof(double));
-    memcpy(vec_tqek1 + offset, rhs.vec_tqek1, chunk_size * sizeof(double));
-    memcpy(vec_tqsk1 + offset, rhs.vec_tqsk1, chunk_size * sizeof(double));
-    memcpy(vec_tqek2 + offset, rhs.vec_tqek2, chunk_size * sizeof(double));
-    memcpy(vec_tqsk2 + offset, rhs.vec_tqsk2, chunk_size * sizeof(double));
-    // Insert vectors of multiple directions per configuration
-    offset = (xi1 - _first_xi) * configurations * ndirs;
-    chunk_size = rhs.xi_block_size * configurations * ndirs;
-    memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_muls + 16 * offset, rhs.vec_dir_muls, 16 * chunk_size * sizeof(double));
-    memcpy(vec_dir_mulslr + 16 * offset, rhs.vec_dir_mulslr, 16 * chunk_size * sizeof(double));
-    // Insert vectors whose sizes depend on wavelengths and directions
-    offset = (xi1 - _first_xi) * ndirs;
-    chunk_size = rhs.xi_block_size * ndirs;
-    memcpy(vec_dir_sat11 + offset, rhs.vec_dir_sat11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sat21 + offset, rhs.vec_dir_sat21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sat12 + offset, rhs.vec_dir_sat12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sat22 + offset, rhs.vec_dir_sat22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_scc1 + offset, rhs.vec_dir_scc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_scc2 + offset, rhs.vec_dir_scc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_abc1 + offset, rhs.vec_dir_abc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_abc2 + offset, rhs.vec_dir_abc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_exc1 + offset, rhs.vec_dir_exc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_exc2 + offset, rhs.vec_dir_exc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_albedc1 + offset, rhs.vec_dir_albedc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_albedc2 + offset, rhs.vec_dir_albedc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_qscc1 + offset, rhs.vec_dir_qscc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_qscc2 + offset, rhs.vec_dir_qscc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_qabc1 + offset, rhs.vec_dir_qabc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_qabc2 + offset, rhs.vec_dir_qabc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_qexc1 + offset, rhs.vec_dir_qexc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_qexc2 + offset, rhs.vec_dir_qexc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_sccrt1 + offset, rhs.vec_dir_sccrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_sccrt2 + offset, rhs.vec_dir_sccrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_abcrt1 + offset, rhs.vec_dir_abcrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_abcrt2 + offset, rhs.vec_dir_abcrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_excrt1 + offset, rhs.vec_dir_excrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_excrt2 + offset, rhs.vec_dir_excrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fsac11 + offset, rhs.vec_dir_fsac11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsac21 + offset, rhs.vec_dir_fsac21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsac12 + offset, rhs.vec_dir_fsac12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsac22 + offset, rhs.vec_dir_fsac22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sac11 + offset, rhs.vec_dir_sac11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sac21 + offset, rhs.vec_dir_sac21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sac12 + offset, rhs.vec_dir_sac12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sac22 + offset, rhs.vec_dir_sac22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_qschuc1 + offset, rhs.vec_dir_qschuc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_qschuc2 + offset, rhs.vec_dir_qschuc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_pschuc1 + offset, rhs.vec_dir_pschuc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_pschuc2 + offset, rhs.vec_dir_pschuc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_s0magc1 + offset, rhs.vec_dir_s0magc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_s0magc2 + offset, rhs.vec_dir_s0magc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_cosavc1 + offset, rhs.vec_dir_cosavc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_cosavc2 + offset, rhs.vec_dir_cosavc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_raprc1 + offset, rhs.vec_dir_raprc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_raprc2 + offset, rhs.vec_dir_raprc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_flc1 + offset, rhs.vec_dir_flc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_flc2 + offset, rhs.vec_dir_flc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_frc1 + offset, rhs.vec_dir_frc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_frc2 + offset, rhs.vec_dir_frc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fkc1 + offset, rhs.vec_dir_fkc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fkc2 + offset, rhs.vec_dir_fkc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fxc1 + offset, rhs.vec_dir_fxc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fxc2 + offset, rhs.vec_dir_fxc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fyc1 + offset, rhs.vec_dir_fyc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fyc2 + offset, rhs.vec_dir_fyc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fzc1 + offset, rhs.vec_dir_fzc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fzc2 + offset, rhs.vec_dir_fzc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqelc1 + offset, rhs.vec_dir_tqelc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqelc2 + offset, rhs.vec_dir_tqelc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqerc1 + offset, rhs.vec_dir_tqerc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqerc2 + offset, rhs.vec_dir_tqerc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqekc1 + offset, rhs.vec_dir_tqekc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqekc2 + offset, rhs.vec_dir_tqekc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqexc1 + offset, rhs.vec_dir_tqexc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqexc2 + offset, rhs.vec_dir_tqexc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqeyc1 + offset, rhs.vec_dir_tqeyc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqeyc2 + offset, rhs.vec_dir_tqeyc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqezc1 + offset, rhs.vec_dir_tqezc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqezc2 + offset, rhs.vec_dir_tqezc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqslc1 + offset, rhs.vec_dir_tqslc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqslc2 + offset, rhs.vec_dir_tqslc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsrc1 + offset, rhs.vec_dir_tqsrc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsrc2 + offset, rhs.vec_dir_tqsrc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqskc1 + offset, rhs.vec_dir_tqskc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqskc2 + offset, rhs.vec_dir_tqskc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsxc1 + offset, rhs.vec_dir_tqsxc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsxc2 + offset, rhs.vec_dir_tqsxc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsyc1 + offset, rhs.vec_dir_tqsyc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsyc2 + offset, rhs.vec_dir_tqsyc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqszc1 + offset, rhs.vec_dir_tqszc1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqszc2 + offset, rhs.vec_dir_tqszc2, chunk_size * sizeof(double));
-    memcpy(vec_dir_mulc + 16 * offset, rhs.vec_dir_mulc, 16 * chunk_size * sizeof(double));
-    memcpy(vec_dir_mulclr + 16 * offset, rhs.vec_dir_mulclr, 16 * chunk_size * sizeof(double));
+  if (rhs.skip_flag != 1) {
+    result += (rhs.nsph == nsph) ? 0 : 1;
+    result += (rhs.inpol == inpol) ? 0 : 1;
+    result += (rhs.iavm == iavm) ? 0 : 1;
+    result += (rhs.isam == isam) ? 0 : 1;
+    result += (rhs._num_theta == _num_theta) ? 0 : 1;
+    result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
+    result += (rhs._num_phi == _num_phi) ? 0 : 1;
+    result += (rhs._num_phis == _num_phis) ? 0 : 1;
+    result += (rhs.ndirs == ndirs) ? 0 : 1;
+    result += (rhs.exri == exri) ? 0 : 1;
+    result += (rhs.idfc == idfc) ? 0 : 1;
+    result += (rhs.configurations == configurations) ? 0 : 1;
+    if (result == 0) {
+      int offset, chunk_size, xi1;
+      xi1 = rhs._first_xi;
+      // Insert vectors whose size depends on wavelengths
+      offset = xi1 - _first_xi;
+      chunk_size = rhs.xi_block_size;
+      memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
+      memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
+      memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
+      memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
+      memcpy(vec_fsat + offset, rhs.vec_fsat, chunk_size * sizeof(dcomplex));
+      memcpy(vec_qschut + offset, rhs.vec_qschut, chunk_size * sizeof(double));
+      memcpy(vec_pschut + offset, rhs.vec_pschut, chunk_size * sizeof(double));
+      memcpy(vec_s0magt + offset, rhs.vec_s0magt, chunk_size * sizeof(double));
+      memcpy(vec_scc1 + offset, rhs.vec_scc1, chunk_size * sizeof(double));
+      memcpy(vec_scc2 + offset, rhs.vec_scc2, chunk_size * sizeof(double));
+      memcpy(vec_abc1 + offset, rhs.vec_abc1, chunk_size * sizeof(double));
+      memcpy(vec_abc2 + offset, rhs.vec_abc2, chunk_size * sizeof(double));
+      memcpy(vec_exc1 + offset, rhs.vec_exc1, chunk_size * sizeof(double));
+      memcpy(vec_exc2 + offset, rhs.vec_exc2, chunk_size * sizeof(double));
+      memcpy(vec_albedc1 + offset, rhs.vec_albedc1, chunk_size * sizeof(double));
+      memcpy(vec_albedc2 + offset, rhs.vec_albedc2, chunk_size * sizeof(double));
+      memcpy(vec_qscamc1 + offset, rhs.vec_qscamc1, chunk_size * sizeof(double));
+      memcpy(vec_qscamc2 + offset, rhs.vec_qscamc2, chunk_size * sizeof(double));
+      memcpy(vec_qabsmc1 + offset, rhs.vec_qabsmc1, chunk_size * sizeof(double));
+      memcpy(vec_qabsmc2 + offset, rhs.vec_qabsmc2, chunk_size * sizeof(double));
+      memcpy(vec_qextmc1 + offset, rhs.vec_qextmc1, chunk_size * sizeof(double));
+      memcpy(vec_qextmc2 + offset, rhs.vec_qextmc2, chunk_size * sizeof(double));
+      memcpy(vec_sccrt1 + offset, rhs.vec_sccrt1, chunk_size * sizeof(double));
+      memcpy(vec_sccrt2 + offset, rhs.vec_sccrt2, chunk_size * sizeof(double));
+      memcpy(vec_abcrt1 + offset, rhs.vec_abcrt1, chunk_size * sizeof(double));
+      memcpy(vec_abcrt2 + offset, rhs.vec_abcrt2, chunk_size * sizeof(double));
+      memcpy(vec_excrt1 + offset, rhs.vec_excrt1, chunk_size * sizeof(double));
+      memcpy(vec_excrt2 + offset, rhs.vec_excrt2, chunk_size * sizeof(double));
+      memcpy(vec_fsac11 + offset, rhs.vec_fsac11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsac21 + offset, rhs.vec_fsac21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsac22 + offset, rhs.vec_fsac22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsac12 + offset, rhs.vec_fsac12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_qschuc1 + offset, rhs.vec_qschuc1, chunk_size * sizeof(double));
+      memcpy(vec_qschuc2 + offset, rhs.vec_qschuc2, chunk_size * sizeof(double));
+      memcpy(vec_pschuc1 + offset, rhs.vec_pschuc1, chunk_size * sizeof(double));
+      memcpy(vec_pschuc2 + offset, rhs.vec_pschuc2, chunk_size * sizeof(double));
+      memcpy(vec_s0magc1 + offset, rhs.vec_s0magc1, chunk_size * sizeof(double));
+      memcpy(vec_s0magc2 + offset, rhs.vec_s0magc2, chunk_size * sizeof(double));
+      memcpy(vec_cosavc1 + offset, rhs.vec_cosavc1, chunk_size * sizeof(double));
+      memcpy(vec_cosavc2 + offset, rhs.vec_cosavc2, chunk_size * sizeof(double));
+      memcpy(vec_raprc1 + offset, rhs.vec_raprc1, chunk_size * sizeof(double));
+      memcpy(vec_raprc2 + offset, rhs.vec_raprc2, chunk_size * sizeof(double));
+      memcpy(vec_fkc1 + offset, rhs.vec_fkc1, chunk_size * sizeof(double));
+      memcpy(vec_fkc2 + offset, rhs.vec_fkc2, chunk_size * sizeof(double));
+      // Insert vectors of multiple configuration values per scale
+      offset = (xi1 - _first_xi) * configurations;
+      chunk_size = rhs.xi_block_size * configurations;
+      memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
+      memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
+      memcpy(vec_sphere_scs + offset, rhs.vec_sphere_scs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_abs + offset, rhs.vec_sphere_abs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_exs + offset, rhs.vec_sphere_exs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_albs + offset, rhs.vec_sphere_albs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_sqscs + offset, rhs.vec_sphere_sqscs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_sqabs + offset, rhs.vec_sphere_sqabs, chunk_size * sizeof(double));
+      memcpy(vec_sphere_sqexs + offset, rhs.vec_sphere_sqexs, chunk_size * sizeof(double));
+      memcpy(vec_fsas + offset, rhs.vec_fsas, chunk_size * sizeof(dcomplex));
+      memcpy(vec_qschus + offset, rhs.vec_qschus, chunk_size * sizeof(double));
+      memcpy(vec_pschus + offset, rhs.vec_pschus, chunk_size * sizeof(double));
+      memcpy(vec_s0mags + offset, rhs.vec_s0mags, chunk_size * sizeof(double));
+      memcpy(vec_cosavs + offset, rhs.vec_cosavs, chunk_size * sizeof(double));
+      memcpy(vec_raprs + offset, rhs.vec_raprs, chunk_size * sizeof(double));
+      memcpy(vec_tqek1 + offset, rhs.vec_tqek1, chunk_size * sizeof(double));
+      memcpy(vec_tqsk1 + offset, rhs.vec_tqsk1, chunk_size * sizeof(double));
+      memcpy(vec_tqek2 + offset, rhs.vec_tqek2, chunk_size * sizeof(double));
+      memcpy(vec_tqsk2 + offset, rhs.vec_tqsk2, chunk_size * sizeof(double));
+      // Insert vectors of multiple directions per configuration
+      offset = (xi1 - _first_xi) * configurations * ndirs;
+      chunk_size = rhs.xi_block_size * configurations * ndirs;
+      memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_muls + 16 * offset, rhs.vec_dir_muls, 16 * chunk_size * sizeof(double));
+      memcpy(vec_dir_mulslr + 16 * offset, rhs.vec_dir_mulslr, 16 * chunk_size * sizeof(double));
+      // Insert vectors whose sizes depend on wavelengths and directions
+      offset = (xi1 - _first_xi) * ndirs;
+      chunk_size = rhs.xi_block_size * ndirs;
+      memcpy(vec_dir_sat11 + offset, rhs.vec_dir_sat11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sat21 + offset, rhs.vec_dir_sat21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sat12 + offset, rhs.vec_dir_sat12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sat22 + offset, rhs.vec_dir_sat22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_scc1 + offset, rhs.vec_dir_scc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_scc2 + offset, rhs.vec_dir_scc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_abc1 + offset, rhs.vec_dir_abc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_abc2 + offset, rhs.vec_dir_abc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_exc1 + offset, rhs.vec_dir_exc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_exc2 + offset, rhs.vec_dir_exc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_albedc1 + offset, rhs.vec_dir_albedc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_albedc2 + offset, rhs.vec_dir_albedc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_qscc1 + offset, rhs.vec_dir_qscc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_qscc2 + offset, rhs.vec_dir_qscc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_qabc1 + offset, rhs.vec_dir_qabc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_qabc2 + offset, rhs.vec_dir_qabc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_qexc1 + offset, rhs.vec_dir_qexc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_qexc2 + offset, rhs.vec_dir_qexc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_sccrt1 + offset, rhs.vec_dir_sccrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_sccrt2 + offset, rhs.vec_dir_sccrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_abcrt1 + offset, rhs.vec_dir_abcrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_abcrt2 + offset, rhs.vec_dir_abcrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_excrt1 + offset, rhs.vec_dir_excrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_excrt2 + offset, rhs.vec_dir_excrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fsac11 + offset, rhs.vec_dir_fsac11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsac21 + offset, rhs.vec_dir_fsac21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsac12 + offset, rhs.vec_dir_fsac12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsac22 + offset, rhs.vec_dir_fsac22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sac11 + offset, rhs.vec_dir_sac11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sac21 + offset, rhs.vec_dir_sac21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sac12 + offset, rhs.vec_dir_sac12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sac22 + offset, rhs.vec_dir_sac22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_qschuc1 + offset, rhs.vec_dir_qschuc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_qschuc2 + offset, rhs.vec_dir_qschuc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_pschuc1 + offset, rhs.vec_dir_pschuc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_pschuc2 + offset, rhs.vec_dir_pschuc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_s0magc1 + offset, rhs.vec_dir_s0magc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_s0magc2 + offset, rhs.vec_dir_s0magc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_cosavc1 + offset, rhs.vec_dir_cosavc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_cosavc2 + offset, rhs.vec_dir_cosavc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_raprc1 + offset, rhs.vec_dir_raprc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_raprc2 + offset, rhs.vec_dir_raprc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_flc1 + offset, rhs.vec_dir_flc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_flc2 + offset, rhs.vec_dir_flc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_frc1 + offset, rhs.vec_dir_frc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_frc2 + offset, rhs.vec_dir_frc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fkc1 + offset, rhs.vec_dir_fkc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fkc2 + offset, rhs.vec_dir_fkc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fxc1 + offset, rhs.vec_dir_fxc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fxc2 + offset, rhs.vec_dir_fxc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fyc1 + offset, rhs.vec_dir_fyc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fyc2 + offset, rhs.vec_dir_fyc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fzc1 + offset, rhs.vec_dir_fzc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fzc2 + offset, rhs.vec_dir_fzc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqelc1 + offset, rhs.vec_dir_tqelc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqelc2 + offset, rhs.vec_dir_tqelc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqerc1 + offset, rhs.vec_dir_tqerc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqerc2 + offset, rhs.vec_dir_tqerc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqekc1 + offset, rhs.vec_dir_tqekc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqekc2 + offset, rhs.vec_dir_tqekc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqexc1 + offset, rhs.vec_dir_tqexc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqexc2 + offset, rhs.vec_dir_tqexc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqeyc1 + offset, rhs.vec_dir_tqeyc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqeyc2 + offset, rhs.vec_dir_tqeyc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqezc1 + offset, rhs.vec_dir_tqezc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqezc2 + offset, rhs.vec_dir_tqezc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqslc1 + offset, rhs.vec_dir_tqslc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqslc2 + offset, rhs.vec_dir_tqslc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsrc1 + offset, rhs.vec_dir_tqsrc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsrc2 + offset, rhs.vec_dir_tqsrc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqskc1 + offset, rhs.vec_dir_tqskc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqskc2 + offset, rhs.vec_dir_tqskc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsxc1 + offset, rhs.vec_dir_tqsxc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsxc2 + offset, rhs.vec_dir_tqsxc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsyc1 + offset, rhs.vec_dir_tqsyc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsyc2 + offset, rhs.vec_dir_tqsyc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqszc1 + offset, rhs.vec_dir_tqszc1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqszc2 + offset, rhs.vec_dir_tqszc2, chunk_size * sizeof(double));
+      memcpy(vec_dir_mulc + 16 * offset, rhs.vec_dir_mulc, 16 * chunk_size * sizeof(double));
+      memcpy(vec_dir_mulclr + 16 * offset, rhs.vec_dir_mulclr, 16 * chunk_size * sizeof(double));
+    }
   }
   return result;
 }
@@ -2356,13 +2374,13 @@ int ClusterOutputInfo::write_legacy(const std::string &output) {
 #ifdef MPI_VERSION
 int ClusterOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
   int result = 0;
-  int skip_flag;
+  int 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(&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(&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 (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);
@@ -2562,191 +2580,200 @@ int ClusterOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
   return result;
 }
 
-
 int ClusterOutputInfo::mpisend(const mixMPI *mpidata) {
   int result = 0;
   int chunk_size;
-  // Send output metadata for configuration cross-check
-  MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  // Wait for process 0 to cross-check the configuration
-  MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  if (result == 0) {
-    // Process 0 confirmed the consistency of configuration. Send the data.
-    // Send vectors of single values per scale
-    MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsat, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0magt, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_albedc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_albedc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qscamc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qscamc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qabsmc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qabsmc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qextmc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qextmc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sccrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sccrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abcrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abcrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_excrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_excrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsac11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsac21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsac22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsac12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschuc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschuc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschuc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschuc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0magc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0magc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosavc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosavc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fkc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fkc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+  if (_skip_flag == 1) {
+    // tell the receiver we are not sending anything
+    int flag = 1;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  }
+  else {
+    // tell the receiver we are sending actual stuff
+    int flag = 0;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Send output metadata for configuration cross-check
+    MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Wait for process 0 to cross-check the configuration
+    MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    if (result == 0) {
+      // Process 0 confirmed the consistency of configuration. Send the data.
+      // Send vectors of single values per scale
+      MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsat, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0magt, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_albedc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_albedc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qscamc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qscamc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qabsmc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qabsmc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qextmc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qextmc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sccrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sccrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abcrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abcrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_excrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_excrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsac11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsac21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsac22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsac12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschuc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschuc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschuc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschuc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0magc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0magc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosavc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosavc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fkc1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fkc2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors of multiple configuration values per scale
-    chunk_size = xi_block_size * configurations;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_scs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_abs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_exs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_albs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sqscs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sqabs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sqexs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschus, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschus, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0mags, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosavs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors of multiple configuration values per scale
+      chunk_size = xi_block_size * configurations;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_scs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_abs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_exs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_albs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sqscs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sqabs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sqexs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschus, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschus, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0mags, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosavs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on directions and configurations.
-    chunk_size = ndirs * configurations * xi_block_size;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_muls, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mulslr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on directions and configurations.
+      chunk_size = ndirs * configurations * xi_block_size;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_muls, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mulslr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on directions and scales.
-    chunk_size = xi_block_size * ndirs;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sat11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sat21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sat12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sat22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_albedc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_albedc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qscc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qscc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qabc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qabc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qexc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qexc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sccrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sccrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abcrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abcrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_excrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_excrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsac11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsac21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsac12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsac22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sac11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sac21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sac12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sac22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qschuc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qschuc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_pschuc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_pschuc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_s0magc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_s0magc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_cosavc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_cosavc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_raprc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_raprc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_flc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_flc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_frc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_frc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fkc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fkc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fxc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fxc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fzc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fzc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqelc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqelc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqerc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqerc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqekc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqekc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqexc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqexc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqeyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqeyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqezc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqezc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqslc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqslc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsrc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsrc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqskc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqskc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsxc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsxc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqszc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqszc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mulc, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mulclr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on directions and scales.
+      chunk_size = xi_block_size * ndirs;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sat11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sat21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sat12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sat22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_albedc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_albedc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qscc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qscc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qabc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qabc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qexc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qexc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sccrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sccrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abcrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abcrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_excrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_excrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsac11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsac21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsac12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsac22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sac11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sac21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sac12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sac22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qschuc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qschuc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_pschuc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_pschuc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_s0magc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_s0magc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_cosavc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_cosavc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_raprc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_raprc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_flc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_flc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_frc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_frc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fkc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fkc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fxc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fxc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fzc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fzc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqelc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqelc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqerc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqerc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqekc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqekc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqexc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqexc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqeyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqeyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqezc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqezc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqslc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqslc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsrc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsrc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqskc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqskc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsxc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsxc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsyc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsyc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqszc1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqszc2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mulc, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mulclr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    }
   }
   return result;
 }
@@ -2758,6 +2785,7 @@ InclusionOutputInfo::InclusionOutputInfo(
   ScattererConfiguration *sc, GeometryConfiguration *gc,
   const mixMPI *mpidata, int first_xi, int xi_length
 ) {
+  _skip_flag = 0;
   nsph = gc->number_of_spheres;
   li = gc->li;
   le = gc->le;
@@ -2946,6 +2974,7 @@ InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_name) {
   HDFFile *hdf_file = new HDFFile(hdf5_name, flags);
   herr_t status = hdf_file->get_status();
   string str_name, str_type;
+  _skip_flag = 0;
   if (status == 0) {
     status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
     status = hdf_file->read("LI", "INT32_(1)", &li);
@@ -3260,127 +3289,141 @@ InclusionOutputInfo::InclusionOutputInfo(const std::string &hdf5_name) {
   }
 }
 
+InclusionOutputInfo::InclusionOutputInfo(const int flag) {
+  /*
+    create a dummy placeholder just to know I should skip MPI_Send and MPI_Recv
+  */
+  if (flag == 1) {
+    _skip_flag = 1;
+  } else {
+    UnrecognizedOutputInfo ex(flag);
+    throw ex;
+  }
+}
+
 InclusionOutputInfo::~InclusionOutputInfo() {
-  delete[] vec_x_coords;
-  delete[] vec_y_coords;
-  delete[] vec_z_coords;
-  delete[] vec_jxi;
-  delete[] vec_ier;
-  delete[] vec_vk;
-  delete[] vec_xi;
-  delete[] vec_sphere_sizes;
-  delete[] vec_sphere_ref_indices;
-  delete[] vec_scs1;
-  delete[] vec_scs2;
-  delete[] vec_abs1;
-  delete[] vec_abs2;
-  delete[] vec_exs1;
-  delete[] vec_exs2;
-  delete[] vec_albeds1;
-  delete[] vec_albeds2;
-  delete[] vec_scsrt1;
-  delete[] vec_scsrt2;
-  delete[] vec_absrt1;
-  delete[] vec_absrt2;
-  delete[] vec_exsrt1;
-  delete[] vec_exsrt2;
-  delete[] vec_qschu1;
-  delete[] vec_qschu2;
-  delete[] vec_pschu1;
-  delete[] vec_pschu2;
-  delete[] vec_s0mag1;
-  delete[] vec_s0mag2;
-  delete[] vec_cosav1;
-  delete[] vec_cosav2;
-  delete[] vec_raprs1;
-  delete[] vec_raprs2;
-  delete[] vec_fk1;
-  delete[] vec_fk2;
-  delete[] vec_fsas11;
-  delete[] vec_fsas21;
-  delete[] vec_fsas22;
-  delete[] vec_fsas12;
-  delete[] vec_dir_tidg;
-  delete[] vec_dir_pidg;
-  delete[] vec_dir_tsdg;
-  delete[] vec_dir_psdg;
-  delete[] vec_dir_scand;
-  delete[] vec_dir_cfmp;
-  delete[] vec_dir_sfmp;
-  delete[] vec_dir_cfsp;
-  delete[] vec_dir_sfsp;
-  delete[] vec_dir_un;
-  delete[] vec_dir_uns;
-  delete[] vec_dir_scs1;
-  delete[] vec_dir_scs2;
-  delete[] vec_dir_abs1;
-  delete[] vec_dir_abs2;
-  delete[] vec_dir_exs1;
-  delete[] vec_dir_exs2;
-  delete[] vec_dir_albeds1;
-  delete[] vec_dir_albeds2;
-  delete[] vec_dir_scsrt1;
-  delete[] vec_dir_scsrt2;
-  delete[] vec_dir_absrt1;
-  delete[] vec_dir_absrt2;
-  delete[] vec_dir_exsrt1;
-  delete[] vec_dir_exsrt2;
-  delete[] vec_dir_fsas11;
-  delete[] vec_dir_fsas21;
-  delete[] vec_dir_fsas12;
-  delete[] vec_dir_fsas22;
-  delete[] vec_dir_sas11;
-  delete[] vec_dir_sas21;
-  delete[] vec_dir_sas12;
-  delete[] vec_dir_sas22;
-  delete[] vec_dir_qschu1;
-  delete[] vec_dir_qschu2;
-  delete[] vec_dir_pschu1;
-  delete[] vec_dir_pschu2;
-  delete[] vec_dir_s0mag1;
-  delete[] vec_dir_s0mag2;
-  delete[] vec_dir_cosav1;
-  delete[] vec_dir_cosav2;
-  delete[] vec_dir_rapr1;
-  delete[] vec_dir_rapr2;
-  delete[] vec_dir_fl1;
-  delete[] vec_dir_fl2;
-  delete[] vec_dir_fr1;
-  delete[] vec_dir_fr2;
-  delete[] vec_dir_fk1;
-  delete[] vec_dir_fk2;
-  delete[] vec_dir_fx1;
-  delete[] vec_dir_fx2;
-  delete[] vec_dir_fy1;
-  delete[] vec_dir_fy2;
-  delete[] vec_dir_fz1;
-  delete[] vec_dir_fz2;
-  delete[] vec_dir_tqel1;
-  delete[] vec_dir_tqel2;
-  delete[] vec_dir_tqer1;
-  delete[] vec_dir_tqer2;
-  delete[] vec_dir_tqek1;
-  delete[] vec_dir_tqek2;
-  delete[] vec_dir_tqex1;
-  delete[] vec_dir_tqex2;
-  delete[] vec_dir_tqey1;
-  delete[] vec_dir_tqey2;
-  delete[] vec_dir_tqez1;
-  delete[] vec_dir_tqez2;
-  delete[] vec_dir_tqsl1;
-  delete[] vec_dir_tqsl2;
-  delete[] vec_dir_tqsr1;
-  delete[] vec_dir_tqsr2;
-  delete[] vec_dir_tqsk1;
-  delete[] vec_dir_tqsk2;
-  delete[] vec_dir_tqsx1;
-  delete[] vec_dir_tqsx2;
-  delete[] vec_dir_tqsy1;
-  delete[] vec_dir_tqsy2;
-  delete[] vec_dir_tqsz1;
-  delete[] vec_dir_tqsz2;
-  delete[] vec_dir_mull;
-  delete[] vec_dir_mulllr;
+  if (_skip_flag != 1) {
+    delete[] vec_x_coords;
+    delete[] vec_y_coords;
+    delete[] vec_z_coords;
+    delete[] vec_jxi;
+    delete[] vec_ier;
+    delete[] vec_vk;
+    delete[] vec_xi;
+    delete[] vec_sphere_sizes;
+    delete[] vec_sphere_ref_indices;
+    delete[] vec_scs1;
+    delete[] vec_scs2;
+    delete[] vec_abs1;
+    delete[] vec_abs2;
+    delete[] vec_exs1;
+    delete[] vec_exs2;
+    delete[] vec_albeds1;
+    delete[] vec_albeds2;
+    delete[] vec_scsrt1;
+    delete[] vec_scsrt2;
+    delete[] vec_absrt1;
+    delete[] vec_absrt2;
+    delete[] vec_exsrt1;
+    delete[] vec_exsrt2;
+    delete[] vec_qschu1;
+    delete[] vec_qschu2;
+    delete[] vec_pschu1;
+    delete[] vec_pschu2;
+    delete[] vec_s0mag1;
+    delete[] vec_s0mag2;
+    delete[] vec_cosav1;
+    delete[] vec_cosav2;
+    delete[] vec_raprs1;
+    delete[] vec_raprs2;
+    delete[] vec_fk1;
+    delete[] vec_fk2;
+    delete[] vec_fsas11;
+    delete[] vec_fsas21;
+    delete[] vec_fsas22;
+    delete[] vec_fsas12;
+    delete[] vec_dir_tidg;
+    delete[] vec_dir_pidg;
+    delete[] vec_dir_tsdg;
+    delete[] vec_dir_psdg;
+    delete[] vec_dir_scand;
+    delete[] vec_dir_cfmp;
+    delete[] vec_dir_sfmp;
+    delete[] vec_dir_cfsp;
+    delete[] vec_dir_sfsp;
+    delete[] vec_dir_un;
+    delete[] vec_dir_uns;
+    delete[] vec_dir_scs1;
+    delete[] vec_dir_scs2;
+    delete[] vec_dir_abs1;
+    delete[] vec_dir_abs2;
+    delete[] vec_dir_exs1;
+    delete[] vec_dir_exs2;
+    delete[] vec_dir_albeds1;
+    delete[] vec_dir_albeds2;
+    delete[] vec_dir_scsrt1;
+    delete[] vec_dir_scsrt2;
+    delete[] vec_dir_absrt1;
+    delete[] vec_dir_absrt2;
+    delete[] vec_dir_exsrt1;
+    delete[] vec_dir_exsrt2;
+    delete[] vec_dir_fsas11;
+    delete[] vec_dir_fsas21;
+    delete[] vec_dir_fsas12;
+    delete[] vec_dir_fsas22;
+    delete[] vec_dir_sas11;
+    delete[] vec_dir_sas21;
+    delete[] vec_dir_sas12;
+    delete[] vec_dir_sas22;
+    delete[] vec_dir_qschu1;
+    delete[] vec_dir_qschu2;
+    delete[] vec_dir_pschu1;
+    delete[] vec_dir_pschu2;
+    delete[] vec_dir_s0mag1;
+    delete[] vec_dir_s0mag2;
+    delete[] vec_dir_cosav1;
+    delete[] vec_dir_cosav2;
+    delete[] vec_dir_rapr1;
+    delete[] vec_dir_rapr2;
+    delete[] vec_dir_fl1;
+    delete[] vec_dir_fl2;
+    delete[] vec_dir_fr1;
+    delete[] vec_dir_fr2;
+    delete[] vec_dir_fk1;
+    delete[] vec_dir_fk2;
+    delete[] vec_dir_fx1;
+    delete[] vec_dir_fx2;
+    delete[] vec_dir_fy1;
+    delete[] vec_dir_fy2;
+    delete[] vec_dir_fz1;
+    delete[] vec_dir_fz2;
+    delete[] vec_dir_tqel1;
+    delete[] vec_dir_tqel2;
+    delete[] vec_dir_tqer1;
+    delete[] vec_dir_tqer2;
+    delete[] vec_dir_tqek1;
+    delete[] vec_dir_tqek2;
+    delete[] vec_dir_tqex1;
+    delete[] vec_dir_tqex2;
+    delete[] vec_dir_tqey1;
+    delete[] vec_dir_tqey2;
+    delete[] vec_dir_tqez1;
+    delete[] vec_dir_tqez2;
+    delete[] vec_dir_tqsl1;
+    delete[] vec_dir_tqsl2;
+    delete[] vec_dir_tqsr1;
+    delete[] vec_dir_tqsr2;
+    delete[] vec_dir_tqsk1;
+    delete[] vec_dir_tqsk2;
+    delete[] vec_dir_tqsx1;
+    delete[] vec_dir_tqsx2;
+    delete[] vec_dir_tqsy1;
+    delete[] vec_dir_tqsy2;
+    delete[] vec_dir_tqsz1;
+    delete[] vec_dir_tqsz2;
+    delete[] vec_dir_mull;
+    delete[] vec_dir_mulllr;
+  }
 }
 
 long InclusionOutputInfo::compute_size() {
@@ -3439,138 +3482,140 @@ long InclusionOutputInfo::compute_size(
 
 int InclusionOutputInfo::insert(const InclusionOutputInfo &rhs) {
   int result = 0;
-  result += (rhs.nsph == nsph) ? 0 : 1;
-  result += (rhs.inpol == inpol) ? 0 : 1;
-  result += (rhs.iavm == iavm) ? 0 : 1;
-  result += (rhs.isam == isam) ? 0 : 1;
-  result += (rhs._num_theta == _num_theta) ? 0 : 1;
-  result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
-  result += (rhs._num_phi == _num_phi) ? 0 : 1;
-  result += (rhs._num_phis == _num_phis) ? 0 : 1;
-  result += (rhs.ndirs == ndirs) ? 0 : 1;
-  result += (rhs.exri == exri) ? 0 : 1;
-  result += (rhs.idfc == idfc) ? 0 : 1;
-  result += (rhs.configurations == configurations) ? 0 : 1;
-  if (result == 0) {
-    int offset, chunk_size, xi1;
-    xi1 = rhs._first_xi;
-    // Insert vectors whose size depends on wavelengths
-    offset = xi1 - _first_xi;
-    chunk_size = rhs.xi_block_size;
-    memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
-    memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
-    memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
-    memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
-    memcpy(vec_scs1 + offset, rhs.vec_scs1, chunk_size * sizeof(double));
-    memcpy(vec_scs2 + offset, rhs.vec_scs2, chunk_size * sizeof(double));
-    memcpy(vec_abs1 + offset, rhs.vec_abs1, chunk_size * sizeof(double));
-    memcpy(vec_abs2 + offset, rhs.vec_abs2, chunk_size * sizeof(double));
-    memcpy(vec_exs1 + offset, rhs.vec_exs1, chunk_size * sizeof(double));
-    memcpy(vec_exs2 + offset, rhs.vec_exs2, chunk_size * sizeof(double));
-    memcpy(vec_albeds1 + offset, rhs.vec_albeds1, chunk_size * sizeof(double));
-    memcpy(vec_albeds2 + offset, rhs.vec_albeds2, chunk_size * sizeof(double));
-    memcpy(vec_scsrt1 + offset, rhs.vec_scsrt1, chunk_size * sizeof(double));
-    memcpy(vec_scsrt2 + offset, rhs.vec_scsrt2, chunk_size * sizeof(double));
-    memcpy(vec_absrt1 + offset, rhs.vec_absrt1, chunk_size * sizeof(double));
-    memcpy(vec_absrt2 + offset, rhs.vec_absrt2, chunk_size * sizeof(double));
-    memcpy(vec_exsrt1 + offset, rhs.vec_exsrt1, chunk_size * sizeof(double));
-    memcpy(vec_exsrt2 + offset, rhs.vec_exsrt2, chunk_size * sizeof(double));
-    memcpy(vec_qschu1 + offset, rhs.vec_qschu1, chunk_size * sizeof(double));
-    memcpy(vec_qschu2 + offset, rhs.vec_qschu2, chunk_size * sizeof(double));
-    memcpy(vec_pschu1 + offset, rhs.vec_pschu1, chunk_size * sizeof(double));
-    memcpy(vec_pschu2 + offset, rhs.vec_pschu2, chunk_size * sizeof(double));
-    memcpy(vec_s0mag1 + offset, rhs.vec_s0mag1, chunk_size * sizeof(double));
-    memcpy(vec_s0mag2 + offset, rhs.vec_s0mag2, chunk_size * sizeof(double));
-    memcpy(vec_cosav1 + offset, rhs.vec_cosav1, chunk_size * sizeof(double));
-    memcpy(vec_cosav2 + offset, rhs.vec_cosav2, chunk_size * sizeof(double));
-    memcpy(vec_raprs1 + offset, rhs.vec_raprs1, chunk_size * sizeof(double));
-    memcpy(vec_raprs2 + offset, rhs.vec_raprs2, chunk_size * sizeof(double));
-    memcpy(vec_fk1 + offset, rhs.vec_fk1, chunk_size * sizeof(double));
-    memcpy(vec_fk2 + offset, rhs.vec_fk2, chunk_size * sizeof(double));
-    memcpy(vec_fsas11 + offset, rhs.vec_fsas11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsas21 + offset, rhs.vec_fsas21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsas22 + offset, rhs.vec_fsas22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_fsas12 + offset, rhs.vec_fsas12, chunk_size * sizeof(dcomplex));
+  if (rhs.skip_flag != 1) {
+    result += (rhs.nsph == nsph) ? 0 : 1;
+    result += (rhs.inpol == inpol) ? 0 : 1;
+    result += (rhs.iavm == iavm) ? 0 : 1;
+    result += (rhs.isam == isam) ? 0 : 1;
+    result += (rhs._num_theta == _num_theta) ? 0 : 1;
+    result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
+    result += (rhs._num_phi == _num_phi) ? 0 : 1;
+    result += (rhs._num_phis == _num_phis) ? 0 : 1;
+    result += (rhs.ndirs == ndirs) ? 0 : 1;
+    result += (rhs.exri == exri) ? 0 : 1;
+    result += (rhs.idfc == idfc) ? 0 : 1;
+    result += (rhs.configurations == configurations) ? 0 : 1;
+    if (result == 0) {
+      int offset, chunk_size, xi1;
+      xi1 = rhs._first_xi;
+      // Insert vectors whose size depends on wavelengths
+      offset = xi1 - _first_xi;
+      chunk_size = rhs.xi_block_size;
+      memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
+      memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
+      memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
+      memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
+      memcpy(vec_scs1 + offset, rhs.vec_scs1, chunk_size * sizeof(double));
+      memcpy(vec_scs2 + offset, rhs.vec_scs2, chunk_size * sizeof(double));
+      memcpy(vec_abs1 + offset, rhs.vec_abs1, chunk_size * sizeof(double));
+      memcpy(vec_abs2 + offset, rhs.vec_abs2, chunk_size * sizeof(double));
+      memcpy(vec_exs1 + offset, rhs.vec_exs1, chunk_size * sizeof(double));
+      memcpy(vec_exs2 + offset, rhs.vec_exs2, chunk_size * sizeof(double));
+      memcpy(vec_albeds1 + offset, rhs.vec_albeds1, chunk_size * sizeof(double));
+      memcpy(vec_albeds2 + offset, rhs.vec_albeds2, chunk_size * sizeof(double));
+      memcpy(vec_scsrt1 + offset, rhs.vec_scsrt1, chunk_size * sizeof(double));
+      memcpy(vec_scsrt2 + offset, rhs.vec_scsrt2, chunk_size * sizeof(double));
+      memcpy(vec_absrt1 + offset, rhs.vec_absrt1, chunk_size * sizeof(double));
+      memcpy(vec_absrt2 + offset, rhs.vec_absrt2, chunk_size * sizeof(double));
+      memcpy(vec_exsrt1 + offset, rhs.vec_exsrt1, chunk_size * sizeof(double));
+      memcpy(vec_exsrt2 + offset, rhs.vec_exsrt2, chunk_size * sizeof(double));
+      memcpy(vec_qschu1 + offset, rhs.vec_qschu1, chunk_size * sizeof(double));
+      memcpy(vec_qschu2 + offset, rhs.vec_qschu2, chunk_size * sizeof(double));
+      memcpy(vec_pschu1 + offset, rhs.vec_pschu1, chunk_size * sizeof(double));
+      memcpy(vec_pschu2 + offset, rhs.vec_pschu2, chunk_size * sizeof(double));
+      memcpy(vec_s0mag1 + offset, rhs.vec_s0mag1, chunk_size * sizeof(double));
+      memcpy(vec_s0mag2 + offset, rhs.vec_s0mag2, chunk_size * sizeof(double));
+      memcpy(vec_cosav1 + offset, rhs.vec_cosav1, chunk_size * sizeof(double));
+      memcpy(vec_cosav2 + offset, rhs.vec_cosav2, chunk_size * sizeof(double));
+      memcpy(vec_raprs1 + offset, rhs.vec_raprs1, chunk_size * sizeof(double));
+      memcpy(vec_raprs2 + offset, rhs.vec_raprs2, chunk_size * sizeof(double));
+      memcpy(vec_fk1 + offset, rhs.vec_fk1, chunk_size * sizeof(double));
+      memcpy(vec_fk2 + offset, rhs.vec_fk2, chunk_size * sizeof(double));
+      memcpy(vec_fsas11 + offset, rhs.vec_fsas11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsas21 + offset, rhs.vec_fsas21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsas22 + offset, rhs.vec_fsas22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_fsas12 + offset, rhs.vec_fsas12, chunk_size * sizeof(dcomplex));
     
-    // Insert vectors whose sizes depend on configurations
-    offset = (xi1 - _first_xi) * (configurations + 1);
-    chunk_size = rhs.xi_block_size * (configurations + 1);
-    memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
-    memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
+      // Insert vectors whose sizes depend on configurations
+      offset = (xi1 - _first_xi) * (configurations + 1);
+      chunk_size = rhs.xi_block_size * (configurations + 1);
+      memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
+      memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
 
-    // Insert vectors whose sizes depend on directions and wavelengths
-    offset = (xi1 - _first_xi) * ndirs;
-    chunk_size = rhs.xi_block_size * ndirs;
-    memcpy(vec_dir_scs1 + offset, rhs.vec_dir_scs1, chunk_size * sizeof(double));
-    memcpy(vec_dir_scs2 + offset, rhs.vec_dir_scs2, chunk_size * sizeof(double));
-    memcpy(vec_dir_abs1 + offset, rhs.vec_dir_abs1, chunk_size * sizeof(double));
-    memcpy(vec_dir_abs2 + offset, rhs.vec_dir_abs2, chunk_size * sizeof(double));
-    memcpy(vec_dir_exs1 + offset, rhs.vec_dir_exs1, chunk_size * sizeof(double));
-    memcpy(vec_dir_exs2 + offset, rhs.vec_dir_exs2, chunk_size * sizeof(double));
-    memcpy(vec_dir_albeds1 + offset, rhs.vec_dir_albeds1, chunk_size * sizeof(double));
-    memcpy(vec_dir_albeds2 + offset, rhs.vec_dir_albeds2, chunk_size * sizeof(double));
-    memcpy(vec_dir_scsrt1 + offset, rhs.vec_dir_scsrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_scsrt2 + offset, rhs.vec_dir_scsrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_absrt1 + offset, rhs.vec_dir_absrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_absrt2 + offset, rhs.vec_dir_absrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_exsrt1 + offset, rhs.vec_dir_exsrt1, chunk_size * sizeof(double));
-    memcpy(vec_dir_exsrt2 + offset, rhs.vec_dir_exsrt2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fsas11 + offset, rhs.vec_dir_fsas11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsas21 + offset, rhs.vec_dir_fsas21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsas12 + offset, rhs.vec_dir_fsas12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_fsas22 + offset, rhs.vec_dir_fsas22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_qschu1 + offset, rhs.vec_dir_qschu1, chunk_size * sizeof(double));
-    memcpy(vec_dir_qschu2 + offset, rhs.vec_dir_qschu2, chunk_size * sizeof(double));
-    memcpy(vec_dir_pschu1 + offset, rhs.vec_dir_pschu1, chunk_size * sizeof(double));
-    memcpy(vec_dir_pschu2 + offset, rhs.vec_dir_pschu2, chunk_size * sizeof(double));
-    memcpy(vec_dir_s0mag1 + offset, rhs.vec_dir_s0mag1, chunk_size * sizeof(double));
-    memcpy(vec_dir_s0mag2 + offset, rhs.vec_dir_s0mag2, chunk_size * sizeof(double));
-    memcpy(vec_dir_cosav1 + offset, rhs.vec_dir_cosav1, chunk_size * sizeof(double));
-    memcpy(vec_dir_cosav2 + offset, rhs.vec_dir_cosav2, chunk_size * sizeof(double));
-    memcpy(vec_dir_rapr1 + offset, rhs.vec_dir_rapr1, chunk_size * sizeof(double));
-    memcpy(vec_dir_rapr2 + offset, rhs.vec_dir_rapr2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fl1 + offset, rhs.vec_dir_fl1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fl2 + offset, rhs.vec_dir_fl2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fr1 + offset, rhs.vec_dir_fr1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fr2 + offset, rhs.vec_dir_fr2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fk1 + offset, rhs.vec_dir_fk1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fk2 + offset, rhs.vec_dir_fk2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fx1 + offset, rhs.vec_dir_fx1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fx2 + offset, rhs.vec_dir_fx2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fy1 + offset, rhs.vec_dir_fy1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fy2 + offset, rhs.vec_dir_fy2, chunk_size * sizeof(double));
-    memcpy(vec_dir_fz1 + offset, rhs.vec_dir_fz1, chunk_size * sizeof(double));
-    memcpy(vec_dir_fz2 + offset, rhs.vec_dir_fz2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqel1 + offset, rhs.vec_dir_tqel1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqel2 + offset, rhs.vec_dir_tqel2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqer1 + offset, rhs.vec_dir_tqer1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqer2 + offset, rhs.vec_dir_tqer2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqek1 + offset, rhs.vec_dir_tqek1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqek2 + offset, rhs.vec_dir_tqek2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqex1 + offset, rhs.vec_dir_tqex1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqex2 + offset, rhs.vec_dir_tqex2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqey1 + offset, rhs.vec_dir_tqey1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqey2 + offset, rhs.vec_dir_tqey2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqez1 + offset, rhs.vec_dir_tqez1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqez2 + offset, rhs.vec_dir_tqez2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsl1 + offset, rhs.vec_dir_tqsl1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsl2 + offset, rhs.vec_dir_tqsl2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsr1 + offset, rhs.vec_dir_tqsr1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsr2 + offset, rhs.vec_dir_tqsr2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsk1 + offset, rhs.vec_dir_tqsk1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsk2 + offset, rhs.vec_dir_tqsk2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsx1 + offset, rhs.vec_dir_tqsx1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsx2 + offset, rhs.vec_dir_tqsx2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsy1 + offset, rhs.vec_dir_tqsy1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsy2 + offset, rhs.vec_dir_tqsy2, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsz1 + offset, rhs.vec_dir_tqsz1, chunk_size * sizeof(double));
-    memcpy(vec_dir_tqsz2 + offset, rhs.vec_dir_tqsz2, chunk_size * sizeof(double));
-    memcpy(vec_dir_mull + 16 * offset, rhs.vec_dir_mull, 16 * chunk_size * sizeof(double));
-    memcpy(vec_dir_mulllr + 16 * offset, rhs.vec_dir_mulllr, 16 * chunk_size * sizeof(double));
+      // Insert vectors whose sizes depend on directions and wavelengths
+      offset = (xi1 - _first_xi) * ndirs;
+      chunk_size = rhs.xi_block_size * ndirs;
+      memcpy(vec_dir_scs1 + offset, rhs.vec_dir_scs1, chunk_size * sizeof(double));
+      memcpy(vec_dir_scs2 + offset, rhs.vec_dir_scs2, chunk_size * sizeof(double));
+      memcpy(vec_dir_abs1 + offset, rhs.vec_dir_abs1, chunk_size * sizeof(double));
+      memcpy(vec_dir_abs2 + offset, rhs.vec_dir_abs2, chunk_size * sizeof(double));
+      memcpy(vec_dir_exs1 + offset, rhs.vec_dir_exs1, chunk_size * sizeof(double));
+      memcpy(vec_dir_exs2 + offset, rhs.vec_dir_exs2, chunk_size * sizeof(double));
+      memcpy(vec_dir_albeds1 + offset, rhs.vec_dir_albeds1, chunk_size * sizeof(double));
+      memcpy(vec_dir_albeds2 + offset, rhs.vec_dir_albeds2, chunk_size * sizeof(double));
+      memcpy(vec_dir_scsrt1 + offset, rhs.vec_dir_scsrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_scsrt2 + offset, rhs.vec_dir_scsrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_absrt1 + offset, rhs.vec_dir_absrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_absrt2 + offset, rhs.vec_dir_absrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_exsrt1 + offset, rhs.vec_dir_exsrt1, chunk_size * sizeof(double));
+      memcpy(vec_dir_exsrt2 + offset, rhs.vec_dir_exsrt2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fsas11 + offset, rhs.vec_dir_fsas11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsas21 + offset, rhs.vec_dir_fsas21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsas12 + offset, rhs.vec_dir_fsas12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_fsas22 + offset, rhs.vec_dir_fsas22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_qschu1 + offset, rhs.vec_dir_qschu1, chunk_size * sizeof(double));
+      memcpy(vec_dir_qschu2 + offset, rhs.vec_dir_qschu2, chunk_size * sizeof(double));
+      memcpy(vec_dir_pschu1 + offset, rhs.vec_dir_pschu1, chunk_size * sizeof(double));
+      memcpy(vec_dir_pschu2 + offset, rhs.vec_dir_pschu2, chunk_size * sizeof(double));
+      memcpy(vec_dir_s0mag1 + offset, rhs.vec_dir_s0mag1, chunk_size * sizeof(double));
+      memcpy(vec_dir_s0mag2 + offset, rhs.vec_dir_s0mag2, chunk_size * sizeof(double));
+      memcpy(vec_dir_cosav1 + offset, rhs.vec_dir_cosav1, chunk_size * sizeof(double));
+      memcpy(vec_dir_cosav2 + offset, rhs.vec_dir_cosav2, chunk_size * sizeof(double));
+      memcpy(vec_dir_rapr1 + offset, rhs.vec_dir_rapr1, chunk_size * sizeof(double));
+      memcpy(vec_dir_rapr2 + offset, rhs.vec_dir_rapr2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fl1 + offset, rhs.vec_dir_fl1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fl2 + offset, rhs.vec_dir_fl2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fr1 + offset, rhs.vec_dir_fr1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fr2 + offset, rhs.vec_dir_fr2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fk1 + offset, rhs.vec_dir_fk1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fk2 + offset, rhs.vec_dir_fk2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fx1 + offset, rhs.vec_dir_fx1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fx2 + offset, rhs.vec_dir_fx2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fy1 + offset, rhs.vec_dir_fy1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fy2 + offset, rhs.vec_dir_fy2, chunk_size * sizeof(double));
+      memcpy(vec_dir_fz1 + offset, rhs.vec_dir_fz1, chunk_size * sizeof(double));
+      memcpy(vec_dir_fz2 + offset, rhs.vec_dir_fz2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqel1 + offset, rhs.vec_dir_tqel1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqel2 + offset, rhs.vec_dir_tqel2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqer1 + offset, rhs.vec_dir_tqer1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqer2 + offset, rhs.vec_dir_tqer2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqek1 + offset, rhs.vec_dir_tqek1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqek2 + offset, rhs.vec_dir_tqek2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqex1 + offset, rhs.vec_dir_tqex1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqex2 + offset, rhs.vec_dir_tqex2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqey1 + offset, rhs.vec_dir_tqey1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqey2 + offset, rhs.vec_dir_tqey2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqez1 + offset, rhs.vec_dir_tqez1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqez2 + offset, rhs.vec_dir_tqez2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsl1 + offset, rhs.vec_dir_tqsl1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsl2 + offset, rhs.vec_dir_tqsl2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsr1 + offset, rhs.vec_dir_tqsr1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsr2 + offset, rhs.vec_dir_tqsr2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsk1 + offset, rhs.vec_dir_tqsk1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsk2 + offset, rhs.vec_dir_tqsk2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsx1 + offset, rhs.vec_dir_tqsx1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsx2 + offset, rhs.vec_dir_tqsx2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsy1 + offset, rhs.vec_dir_tqsy1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsy2 + offset, rhs.vec_dir_tqsy2, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsz1 + offset, rhs.vec_dir_tqsz1, chunk_size * sizeof(double));
+      memcpy(vec_dir_tqsz2 + offset, rhs.vec_dir_tqsz2, chunk_size * sizeof(double));
+      memcpy(vec_dir_mull + 16 * offset, rhs.vec_dir_mull, 16 * chunk_size * sizeof(double));
+      memcpy(vec_dir_mulllr + 16 * offset, rhs.vec_dir_mulllr, 16 * chunk_size * sizeof(double));
+    }
   }
   return result;
 }
@@ -4485,13 +4530,13 @@ 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 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(&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(&flag, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  // Proceed with the rest _only if_ flag==0, else nothing is to be received
+  if (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);
@@ -4648,140 +4693,150 @@ int InclusionOutputInfo::mpireceive(const mixMPI* mpidata, int pid) {
 int InclusionOutputInfo::mpisend(const mixMPI *mpidata) {
   int result = 0;
   int chunk_size;
-  // Send output metadata for configuration cross-check
-  MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  // Wait for process 0 to cross-check the configuration
-  MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  if (result == 0) {
-    // Process 0 confirmed the consistency of configuration. Send the data.
-    // Send vectors of single values per scale
-    MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_albeds1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_albeds2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_absrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_absrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0mag1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0mag2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosav1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosav2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fk1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fk2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+  if (_skip_flag == 1) {
+    // tell the receiver we are not sending anything
+    int flag = 1;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  }
+  else {
+    // tell the receiver we are sending actual stuff
+    int flag = 0;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Send output metadata for configuration cross-check
+    MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&iavm, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Wait for process 0 to cross-check the configuration
+    MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    if (result == 0) {
+      // Process 0 confirmed the consistency of configuration. Send the data.
+      // Send vectors of single values per scale
+      MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_albeds1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_albeds2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_absrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_absrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exsrt1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exsrt2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschu1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschu2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0mag1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0mag2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosav1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosav2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprs1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprs2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fk1, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fk2, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas11, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas21, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas22, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas12, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on configurations
-    chunk_size = xi_block_size * (configurations + 1);
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on configurations
+      chunk_size = xi_block_size * (configurations + 1);
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on directions and wavelengths
-    chunk_size = xi_block_size * ndirs;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_abs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_albeds1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_albeds2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_scsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_absrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_absrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_exsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fsas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_qschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_pschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_pschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_s0mag1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_s0mag2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_cosav1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_cosav2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_rapr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_rapr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqel1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqel2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqer1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqer2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqex1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqex2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqey1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqey2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqez1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqez2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_tqsz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mull, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mulllr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on directions and wavelengths
+      chunk_size = xi_block_size * ndirs;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_abs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exs1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exs2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_albeds1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_albeds2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_scsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_absrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_absrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exsrt1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_exsrt2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fsas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_qschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_pschu1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_pschu2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_s0mag1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_s0mag2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_cosav1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_cosav2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_rapr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_rapr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqel1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqel2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqer1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqer2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqex1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqex2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqey1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqey2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqez1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqez2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsl1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsl2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsr1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsr2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsx1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsx2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsy1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsy2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsz1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_tqsz2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mull, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mulllr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    }
   }
   return result;
 }
@@ -4793,6 +4848,7 @@ SphereOutputInfo::SphereOutputInfo(
   ScattererConfiguration *sc, GeometryConfiguration *gc,
   const mixMPI *mpidata, int first_xi, int xi_length
 ) {
+  _skip_flag = 0;
   _first_xi = first_xi;
   nsph = gc->number_of_spheres;
   lm = gc->l_max;
@@ -4905,6 +4961,7 @@ SphereOutputInfo::SphereOutputInfo(const std::string &hdf5_name) {
   HDFFile *hdf_file = new HDFFile(hdf5_name, flags);
   herr_t status = hdf_file->get_status();
   string str_name, str_type;
+  _skip_flag = 0;
   if (status == 0) {
     status = hdf_file->read("NSPH", "INT32_(1)", &nsph);
     status = hdf_file->read("LM", "INT32_(1)", &lm);
@@ -5077,56 +5134,70 @@ SphereOutputInfo::SphereOutputInfo(const std::string &hdf5_name) {
   }
 }
 
+SphereOutputInfo::SphereOutputInfo(const int flag) {
+  /*
+    create a dummy placeholder just to know I should skip MPI_Send and MPI_Recv
+  */
+  if (flag == 1) {
+    _skip_flag = 1;
+  } else {
+    UnrecognizedOutputInfo ex(flag);
+    throw ex;
+  }
+}
+
 SphereOutputInfo::~SphereOutputInfo() {
-  delete[] vec_jxi;
-  delete[] vec_ier;
-  delete[] vec_vk;
-  delete[] vec_xi;
-  delete[] vec_sphere_sizes;
-  delete[] vec_sphere_ref_indices;
-  delete[] vec_scs;
-  delete[] vec_abs;
-  delete[] vec_exs;
-  delete[] vec_albeds;
-  delete[] vec_scsrt;
-  delete[] vec_absrt;
-  delete[] vec_exsrt;
-  delete[] vec_fsas;
-  delete[] vec_qschu;
-  delete[] vec_pschu;
-  delete[] vec_s0mag;
-  delete[] vec_cosav;
-  delete[] vec_raprs;
-  delete[] vec_tqek1;
-  delete[] vec_tqek2;
-  delete[] vec_tqsk1;
-  delete[] vec_tqsk2;
-  if (nsph != 1) {
-    delete[] vec_fsat;
-    delete[] vec_qschut;
-    delete[] vec_pschut;
-    delete[] vec_s0magt;
+  if (_skip_flag != 1) {
+    delete[] vec_jxi;
+    delete[] vec_ier;
+    delete[] vec_vk;
+    delete[] vec_xi;
+    delete[] vec_sphere_sizes;
+    delete[] vec_sphere_ref_indices;
+    delete[] vec_scs;
+    delete[] vec_abs;
+    delete[] vec_exs;
+    delete[] vec_albeds;
+    delete[] vec_scsrt;
+    delete[] vec_absrt;
+    delete[] vec_exsrt;
+    delete[] vec_fsas;
+    delete[] vec_qschu;
+    delete[] vec_pschu;
+    delete[] vec_s0mag;
+    delete[] vec_cosav;
+    delete[] vec_raprs;
+    delete[] vec_tqek1;
+    delete[] vec_tqek2;
+    delete[] vec_tqsk1;
+    delete[] vec_tqsk2;
+    if (nsph != 1) {
+      delete[] vec_fsat;
+      delete[] vec_qschut;
+      delete[] vec_pschut;
+      delete[] vec_s0magt;
+    }
+    delete[] vec_dir_tidg;
+    delete[] vec_dir_pidg;
+    delete[] vec_dir_tsdg;
+    delete[] vec_dir_psdg;
+    delete[] vec_dir_scand;
+    delete[] vec_dir_cfmp;
+    delete[] vec_dir_cfsp;
+    delete[] vec_dir_sfmp;
+    delete[] vec_dir_sfsp;
+    delete[] vec_dir_un;
+    delete[] vec_dir_uns;
+    delete[] vec_dir_sas11;
+    delete[] vec_dir_sas21;
+    delete[] vec_dir_sas12;
+    delete[] vec_dir_sas22;
+    delete[] vec_dir_fx;
+    delete[] vec_dir_fy;
+    delete[] vec_dir_fz;
+    delete[] vec_dir_muls;
+    delete[] vec_dir_mulslr;
   }
-  delete[] vec_dir_tidg;
-  delete[] vec_dir_pidg;
-  delete[] vec_dir_tsdg;
-  delete[] vec_dir_psdg;
-  delete[] vec_dir_scand;
-  delete[] vec_dir_cfmp;
-  delete[] vec_dir_cfsp;
-  delete[] vec_dir_sfmp;
-  delete[] vec_dir_sfsp;
-  delete[] vec_dir_un;
-  delete[] vec_dir_uns;
-  delete[] vec_dir_sas11;
-  delete[] vec_dir_sas21;
-  delete[] vec_dir_sas12;
-  delete[] vec_dir_sas22;
-  delete[] vec_dir_fx;
-  delete[] vec_dir_fy;
-  delete[] vec_dir_fz;
-  delete[] vec_dir_muls;
-  delete[] vec_dir_mulslr;
 }
 
 long SphereOutputInfo::compute_size(
@@ -5202,74 +5273,76 @@ long SphereOutputInfo::compute_size() {
 
 int SphereOutputInfo::insert(const SphereOutputInfo &rhs) {
   int result = 0;
-  result += (rhs.nsph == nsph) ? 0 : 1;
-  result += (rhs.inpol == inpol) ? 0 : 1;
-  result += (rhs.isam == isam) ? 0 : 1;
-  result += (rhs._num_theta == _num_theta) ? 0 : 1;
-  result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
-  result += (rhs._num_phi == _num_phi) ? 0 : 1;
-  result += (rhs._num_phis == _num_phis) ? 0 : 1;
-  result += (rhs.ndirs == ndirs) ? 0 : 1;
-  result += (rhs.exri == exri) ? 0 : 1;
-  result += (rhs.idfc == idfc) ? 0 : 1;
-  result += (rhs.configurations == configurations) ? 0 : 1;
-  if (result == 0) {
-    int offset, chunk_size, xi1;
-    xi1 = rhs._first_xi;
-    // Insert vectors whose sizes depend on wavelengths
-    offset = xi1 - _first_xi;
-    chunk_size = rhs.xi_block_size;
-    memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
-    memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
-    memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
-    memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
-    if (nsph != 1) {
-      memcpy(vec_fsat + offset, rhs.vec_fsat, chunk_size * sizeof(dcomplex));
-      memcpy(vec_qschut + offset, rhs.vec_qschut, chunk_size * sizeof(double));
-      memcpy(vec_pschut + offset, rhs.vec_pschut, chunk_size * sizeof(double));
-      memcpy(vec_s0magt + offset, rhs.vec_s0magt, chunk_size * sizeof(double));
-    }
+  if (rhs.skip_flag != 1) {
+    result += (rhs.nsph == nsph) ? 0 : 1;
+    result += (rhs.inpol == inpol) ? 0 : 1;
+    result += (rhs.isam == isam) ? 0 : 1;
+    result += (rhs._num_theta == _num_theta) ? 0 : 1;
+    result += (rhs._num_thetas == _num_thetas) ? 0 : 1;
+    result += (rhs._num_phi == _num_phi) ? 0 : 1;
+    result += (rhs._num_phis == _num_phis) ? 0 : 1;
+    result += (rhs.ndirs == ndirs) ? 0 : 1;
+    result += (rhs.exri == exri) ? 0 : 1;
+    result += (rhs.idfc == idfc) ? 0 : 1;
+    result += (rhs.configurations == configurations) ? 0 : 1;
+    if (result == 0) {
+      int offset, chunk_size, xi1;
+      xi1 = rhs._first_xi;
+      // Insert vectors whose sizes depend on wavelengths
+      offset = xi1 - _first_xi;
+      chunk_size = rhs.xi_block_size;
+      memcpy(vec_jxi + offset, rhs.vec_jxi, chunk_size * sizeof(int));
+      memcpy(vec_ier + offset, rhs.vec_ier, chunk_size * sizeof(short));
+      memcpy(vec_vk + offset, rhs.vec_vk, chunk_size * sizeof(double));
+      memcpy(vec_xi + offset, rhs.vec_xi, chunk_size * sizeof(double));
+      if (nsph != 1) {
+	memcpy(vec_fsat + offset, rhs.vec_fsat, chunk_size * sizeof(dcomplex));
+	memcpy(vec_qschut + offset, rhs.vec_qschut, chunk_size * sizeof(double));
+	memcpy(vec_pschut + offset, rhs.vec_pschut, chunk_size * sizeof(double));
+	memcpy(vec_s0magt + offset, rhs.vec_s0magt, chunk_size * sizeof(double));
+      }
 
-    // Insert vectors whose sizes depend on configurations and wavelengths
-    offset = (xi1 - _first_xi) * configurations;
-    chunk_size = rhs.xi_block_size * configurations;
-    memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
-    memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
-    memcpy(vec_scs + offset, rhs.vec_scs, chunk_size * sizeof(double));
-    memcpy(vec_abs + offset, rhs.vec_abs, chunk_size * sizeof(double));
-    memcpy(vec_exs + offset, rhs.vec_exs, chunk_size * sizeof(double));
-    memcpy(vec_albeds + offset, rhs.vec_albeds, chunk_size * sizeof(double));
-    memcpy(vec_scsrt + offset, rhs.vec_scsrt, chunk_size * sizeof(double));
-    memcpy(vec_absrt + offset, rhs.vec_absrt, chunk_size * sizeof(double));
-    memcpy(vec_exsrt + offset, rhs.vec_exsrt, chunk_size * sizeof(double));
-    memcpy(vec_fsas + offset, rhs.vec_fsas, chunk_size * sizeof(dcomplex));
-    memcpy(vec_qschu + offset, rhs.vec_qschu, chunk_size * sizeof(double));
-    memcpy(vec_pschu + offset, rhs.vec_pschu, chunk_size * sizeof(double));
-    memcpy(vec_s0mag + offset, rhs.vec_s0mag, chunk_size * sizeof(double));
-    memcpy(vec_cosav + offset, rhs.vec_cosav, chunk_size * sizeof(double));
-    memcpy(vec_raprs + offset, rhs.vec_raprs, chunk_size * sizeof(double));
-    memcpy(vec_tqek1 + offset, rhs.vec_tqek1, chunk_size * sizeof(double));
-    memcpy(vec_tqek2 + offset, rhs.vec_tqek2, chunk_size * sizeof(double));
-    memcpy(vec_tqsk1 + offset, rhs.vec_tqsk1, chunk_size * sizeof(double));
-    memcpy(vec_tqsk2 + offset, rhs.vec_tqsk2, chunk_size * sizeof(double));
+      // Insert vectors whose sizes depend on configurations and wavelengths
+      offset = (xi1 - _first_xi) * configurations;
+      chunk_size = rhs.xi_block_size * configurations;
+      memcpy(vec_sphere_sizes + offset, rhs.vec_sphere_sizes, chunk_size * sizeof(double));
+      memcpy(vec_sphere_ref_indices + offset, rhs.vec_sphere_ref_indices, chunk_size * sizeof(dcomplex));
+      memcpy(vec_scs + offset, rhs.vec_scs, chunk_size * sizeof(double));
+      memcpy(vec_abs + offset, rhs.vec_abs, chunk_size * sizeof(double));
+      memcpy(vec_exs + offset, rhs.vec_exs, chunk_size * sizeof(double));
+      memcpy(vec_albeds + offset, rhs.vec_albeds, chunk_size * sizeof(double));
+      memcpy(vec_scsrt + offset, rhs.vec_scsrt, chunk_size * sizeof(double));
+      memcpy(vec_absrt + offset, rhs.vec_absrt, chunk_size * sizeof(double));
+      memcpy(vec_exsrt + offset, rhs.vec_exsrt, chunk_size * sizeof(double));
+      memcpy(vec_fsas + offset, rhs.vec_fsas, chunk_size * sizeof(dcomplex));
+      memcpy(vec_qschu + offset, rhs.vec_qschu, chunk_size * sizeof(double));
+      memcpy(vec_pschu + offset, rhs.vec_pschu, chunk_size * sizeof(double));
+      memcpy(vec_s0mag + offset, rhs.vec_s0mag, chunk_size * sizeof(double));
+      memcpy(vec_cosav + offset, rhs.vec_cosav, chunk_size * sizeof(double));
+      memcpy(vec_raprs + offset, rhs.vec_raprs, chunk_size * sizeof(double));
+      memcpy(vec_tqek1 + offset, rhs.vec_tqek1, chunk_size * sizeof(double));
+      memcpy(vec_tqek2 + offset, rhs.vec_tqek2, chunk_size * sizeof(double));
+      memcpy(vec_tqsk1 + offset, rhs.vec_tqsk1, chunk_size * sizeof(double));
+      memcpy(vec_tqsk2 + offset, rhs.vec_tqsk2, chunk_size * sizeof(double));
     
-    // Insert vectors whose sizes depend on NSPH, directions and wavelengths
-    offset = (xi1 - _first_xi) * nsph * ndirs;
-    chunk_size = rhs.xi_block_size * nsph * ndirs;
-    memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
-    memcpy(vec_dir_muls + 16 * offset, rhs.vec_dir_muls, 16 * chunk_size * sizeof(double));
-    memcpy(vec_dir_mulslr + 16 * offset, rhs.vec_dir_mulslr, 16 * chunk_size * sizeof(double));
+      // Insert vectors whose sizes depend on NSPH, directions and wavelengths
+      offset = (xi1 - _first_xi) * nsph * ndirs;
+      chunk_size = rhs.xi_block_size * nsph * ndirs;
+      memcpy(vec_dir_sas11 + offset, rhs.vec_dir_sas11, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas21 + offset, rhs.vec_dir_sas21, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas12 + offset, rhs.vec_dir_sas12, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_sas22 + offset, rhs.vec_dir_sas22, chunk_size * sizeof(dcomplex));
+      memcpy(vec_dir_muls + 16 * offset, rhs.vec_dir_muls, 16 * chunk_size * sizeof(double));
+      memcpy(vec_dir_mulslr + 16 * offset, rhs.vec_dir_mulslr, 16 * chunk_size * sizeof(double));
 
-    // Insert vectors whose sizes depend on NSPH, incidence directions and wavelengths
-    offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
-    chunk_size = rhs.xi_block_size  * nsph * _num_theta * _num_phi;
-    memcpy(vec_dir_fx + offset, rhs.vec_dir_fx, chunk_size * sizeof(double));
-    memcpy(vec_dir_fy + offset, rhs.vec_dir_fy, chunk_size * sizeof(double));
-    memcpy(vec_dir_fz + offset, rhs.vec_dir_fz, chunk_size * sizeof(double));
-    // TODO: fix the vector sizes in HDF5 writer and MPI communicators
+      // Insert vectors whose sizes depend on NSPH, incidence directions and wavelengths
+      offset = (xi1 - _first_xi) * nsph * _num_theta * _num_phi;
+      chunk_size = rhs.xi_block_size  * nsph * _num_theta * _num_phi;
+      memcpy(vec_dir_fx + offset, rhs.vec_dir_fx, chunk_size * sizeof(double));
+      memcpy(vec_dir_fy + offset, rhs.vec_dir_fy, chunk_size * sizeof(double));
+      memcpy(vec_dir_fz + offset, rhs.vec_dir_fz, chunk_size * sizeof(double));
+      // TODO: fix the vector sizes in HDF5 writer and MPI communicators
+    }
   }
   return result;
 }
@@ -5746,13 +5819,13 @@ 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 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(&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(&flag, 1, MPI_INT32_T, pid, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  // Proceed with the rest _only if_ flag==0, else nothing is to be received
+  if (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);
@@ -5842,74 +5915,84 @@ int SphereOutputInfo::mpireceive(const mixMPI *mpidata, int pid) {
 int SphereOutputInfo::mpisend(const mixMPI *mpidata) {
   int result = 0;
   int chunk_size;
-  // Send output metadata for configuration cross-check
-  MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-  // Wait for process 0 to cross-check the configuration
-  MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-  if (result == 0) {
-    // Process 0 confirmed the consistency of configuration. Send the data.
-    // Send vectors of single values per scale
-    MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    if (nsph != 1) {
-      MPI_Send(vec_fsat, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-      MPI_Send(vec_qschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-      MPI_Send(vec_pschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-      MPI_Send(vec_s0magt, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    }
+  if (_skip_flag == 1) {
+    // tell the receiver we are not sending anything
+    int flag = 1;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+  }
+  else {
+    // tell the receiver we are sending actual stuff
+    int flag = 0;
+    MPI_Send(&flag, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Send output metadata for configuration cross-check
+    MPI_Send(&nsph, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&inpol, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&isam, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_theta, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_thetas, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&_num_phis, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&ndirs, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&exri, 1, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&idfc, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    MPI_Send(&configurations, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+    // Wait for process 0 to cross-check the configuration
+    MPI_Recv(&result, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    if (result == 0) {
+      // Process 0 confirmed the consistency of configuration. Send the data.
+      // Send vectors of single values per scale
+      MPI_Send(&_first_xi, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(&xi_block_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_jxi, xi_block_size, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_ier, xi_block_size, MPI_SHORT, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_vk, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_xi, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      if (nsph != 1) {
+	MPI_Send(vec_fsat, xi_block_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+	MPI_Send(vec_qschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+	MPI_Send(vec_pschut, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+	MPI_Send(vec_s0magt, xi_block_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      }
 
-    // Send vectors whose sizes depend on configurations and scales
-    chunk_size = xi_block_size * configurations;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_abs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_albeds, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_scsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_absrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_exsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_fsas, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_qschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_pschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_s0mag, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_cosav, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_raprs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on configurations and scales
+      chunk_size = xi_block_size * configurations;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_sizes, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_sphere_ref_indices, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_abs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_albeds, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_scsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_absrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_exsrt, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_fsas, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_qschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_pschu, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_s0mag, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_cosav, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_raprs, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqek1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqek2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqsk1, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_tqsk2, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on NSPH, directions and wavelengths
-    chunk_size = xi_block_size * nsph * ndirs;
-    MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_muls, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_mulslr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on NSPH, directions and wavelengths
+      chunk_size = xi_block_size * nsph * ndirs;
+      MPI_Send(&chunk_size, 1, MPI_INT32_T, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas11, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas21, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas12, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_sas22, chunk_size, MPI_C_DOUBLE_COMPLEX, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_muls, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_mulslr, 16 * chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
 
-    // Send vectors whose sizes depend on NSPH, incidence directions and wavelengths
-    chunk_size = xi_block_size * nsph * _num_theta * _num_phi;
-    MPI_Send(vec_dir_fx, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fy, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
-    MPI_Send(vec_dir_fz, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      // Send vectors whose sizes depend on NSPH, incidence directions and wavelengths
+      chunk_size = xi_block_size * nsph * _num_theta * _num_phi;
+      MPI_Send(vec_dir_fx, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fy, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+      MPI_Send(vec_dir_fz, chunk_size, MPI_DOUBLE, 0, 10, MPI_COMM_WORLD);
+    }
   }
   return result;
 }
diff --git a/src/sphere/sphere.cpp b/src/sphere/sphere.cpp
index 18d19b80bc716356e40f39449d21bde53c2af666..288c259f4aee6b48efe3b6980af1e87f01e6dbd8 100644
--- a/src/sphere/sphere.cpp
+++ b/src/sphere/sphere.cpp
@@ -330,19 +330,17 @@ void sphere(const string& config_file, const string& data_file, const string& ou
 	    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 there is no input for this thread, set output pointer to NULL.
-	      p_outarray[myompthread] = NULL;
+	      // If there is no input for this thread, mark to skip.
+	      p_outarray[myompthread] = new SphereOutputInfo(1);
 	    }
 	  }
 #pragma omp barrier
 	  // threads different from 0 append their virtual files to the one of thread 0, and delete them
 	  if (myompthread == 0) {
 	    for (int ti=1; ti<ompnumthreads; ti++) {
-	      if (p_outarray[ti] != NULL) {
-		p_outarray[0]->insert(*(p_outarray[ti]));
-		delete p_outarray[ti];
-		p_outarray[ti] = NULL;
-	      }
+	      p_outarray[0]->insert(*(p_outarray[ti]));
+	      delete p_outarray[ti];
+	      p_outarray[ti] = NULL;
 	      vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	      delete vtppoanarray[ti];
 	    }
@@ -492,39 +490,25 @@ void sphere(const string& config_file, const string& data_file, const string& ou
 	  }
 	  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 there is no input for this thread, set the output pointer to NULL.
-	    p_outarray[myompthread] = NULL;
-	    //}	  
+	  p_outarray[myompthread] = new SphereOutputInfo(1);
 	}
 
 #pragma omp barrier
 	// threads different from 0 append their virtual files to the one of thread 0, and delete them
 	if (myompthread == 0) {
 	  for (int ti=1; ti<ompnumthreads; ti++) {
-	    if (p_outarray[ti] != NULL) {
-	      p_outarray[0]->insert(*(p_outarray[ti]));
-	      delete p_outarray[ti];
-	      p_outarray[ti] = NULL;
-	    }
+	    p_outarray[0]->insert(*(p_outarray[ti]));
+	    delete p_outarray[ti];
+	    p_outarray[ti] = NULL;
 	    vtppoanarray[0]->append(*(vtppoanarray[ti]));
 	    delete vtppoanarray[ti];
 	  }
 	  // 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) {
-	      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];
-	      }
+	      p_outarray[0]->mpisend(mpidata);
+	      delete p_outarray[0];
+	      p_outarray[0] = NULL;
 	      vtppoanarray[0]->mpisend(mpidata);
 	      delete vtppoanarray[0];
 	    }