#pragma once

// [io-base3d.h]: I/O: 3D array
//
// -------------------------------------------------------------------------------------------- //

#define _CRT_SECURE_NO_DEPRECATE
#include <stdio.h>

#include "nse-alloc.h"
#include "grid-id.h"
#include "io-misc.h"

#include "bin-stamp.h"
#include "bin-named-stamp.h"

namespace nse
{
	// Tecplot //
	template< typename T >
	int write_tecplot(const std::string& filename,
		const T* out, const T* cx, const T* cy, const T* cz, // [array, coordinates]
		const int nx, const int ny, const int nz,
		const int ib, const int ie,
		const int jb, const int je,
		const int kb, const int ke,

		const char* title, const char* name,
		const char* name_dimx, const char* name_dimy, const char* name_dimz,
		const T time);

	template< typename T >
	int write_tecplot(const std::string& filename,
		const T* uout, const T* vout, const T* wout,
		const T* cx, const T* cy, const T* cz, // [array, coordinates]
		const int nx, const int ny, const int nz,
		const int ib, const int ie,
		const int jb, const int je,
		const int kb, const int ke,

		const char* title, const char* uname, const char* vname, const char* wname,
		const char* name_dimx, const char* name_dimy, const char* name_dimz,
		const T time);
	// -------------------------------------------------------------------------------------------- //


	// Binary //
	template< typename T >
	int write_binary_stamp(const std::string& filename,
		const binStamp< int >& index_stamp,
		const binStamp< double >& cpu_stamp,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId<T>& id, const T time);

	template< typename T >
	int write_binary_stamp(const std::string& filename,
		const binNamedStamp< int >& index_stamp,
		const binNamedStamp< double >& cpu_stamp,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId<T>& id, const T time);

	template< typename T >
	int write_binary(const std::string& filename,
		const T* xin, const char* name,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId<T>& id, const T time);

	template< typename T >
	int write_binary(const std::string& filename,
		const T* uin, const T* vin, const T* win,
		const char* uname, const char* vname, const char* wname,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId<T>& id, const T time);


	template< typename T >
	int read_binary_stamp(const std::string& filename,
		binStamp< int >& index_stamp,
		binStamp< double >& cpu_stamp,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId<T>& id, T* time);

	template< typename T >
	int read_binary_stamp(const std::string& filename,
		binNamedStamp< int >& index_stamp,
		binNamedStamp< double >& cpu_stamp,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId<T>& id, T* time);

	template< typename T >
	int read_binary(const std::string& filename,
		T** xout, char** name,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId<T>& id, T* time);

	template< typename T >
	int read_binary(const std::string& filename,
		T** uout, T** vout, T** wout,
		char** uname, char** vname, char** wname,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId<T>& id, T* time);
	// -------------------------------------------------------------------------------------------- //

	// MPI-I/O 3D datatype //
	template< typename T >
	void mpi_io_write_datatype(MPI_Datatype* file_view, MPI_Datatype* local_view,

		const int mpi_nx, const int mpi_ny, const int mpi_nz,
		const int nx, const int ny, const int nz,
		const int gcx, const int gcy, const int gcz,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z);

	template< typename T >
	void mpi_io_read_datatype(MPI_Datatype* file_view, MPI_Datatype* local_view,

		const int mpi_nx, const int mpi_ny, const int mpi_nz,
		const int nx, const int ny, const int nz,
		const int gcx, const int gcy, const int gcz,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z);
	// -------------------------------------------------------------------------------------------- //

	// MPI-Binary //
	// [name, coordinates, time] - local, on head-rank
	// [array] - distributed
	// [id] - global, on all ranks in comm
	template< typename T >
	int mpi_write_binary(const std::string& filename,
		const char* mpi_datarep,
		const MPI_Comm comm, const int header_rank,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

		const T* xin, const char* name,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId< T >& id, const T time);

	template< typename T >
	int mpi_write_binary(const std::string& filename,
		const char* mpi_datarep,
		const MPI_Comm comm, const int header_rank,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

		const T* uin, const T* vin, const T* win,
		const char* uname, const char* vname, const char* wname,

		const T* cx, const T* cy, const T* cz,
		const T* ex, const T* ey, const T* ez,
		const GridId< T >& id, const T time);


	template< typename T >
	int mpi_read_binary(const std::string& filename,
		const char* mpi_datarep,
		const MPI_Comm comm, const int header_rank,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

		T** xout, char** name,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId< T >& id, T* time);

	template< typename T >
	int mpi_read_binary(const std::string& filename,
		const char* mpi_datarep,
		const MPI_Comm comm, const int header_rank,
		const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

		T** uout, T** vout, T** wout,
		char** uname, char** vname, char** wname,

		T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
		GridId< T >& id, T* time);
	// -------------------------------------------------------------------------------------------- //
}

// Implementation
// -------------------------------------------------------------------------------------------- //

// Tecplot //
template< typename T >
int nse::write_tecplot(const std::string& filename,
	const T* out, const T* cx, const T* cy, const T* cz, // [array, coordinates]
	const int nx, const int ny, const int nz,
	const int ib, const int ie,
	const int jb, const int je,
	const int kb, const int ke,

	const char* title, const char* name,
	const char* name_dimx, const char* name_dimy, const char* name_dimz,
	const T time)
{
	FILE* ptr = fopen(filename.c_str(), "w");
	if (ptr == NULL) return 0;

	fprintf(ptr, " TITLE = \"%s [F(%s,%s,%s)]\"\n",
		title, name_dimx, name_dimy, name_dimz);
	fprintf(ptr, " VARIABLES = \"%s\", \"%s\", \"%s\", \"%s\"\n",
		name_dimx, name_dimy, name_dimz, name);
	fprintf(ptr, " ZONE I = %i, J = %i, K = %i, DATAPACKING = POINT, SOLUTIONTIME = %f\n",
		ie - ib + 1, je - jb + 1, ke - kb + 1, time);

	const int nyz = ny * nz;
	int i, j, k, idx;

	for (k = kb; k <= ke; k++)
		for (j = jb; j <= je; j++)
			for (i = ib; i <= ie; i++)
			{
				idx = i * nyz + j * nz + k;
				fprintf(ptr, "%f %f %f %f\n", cx[i], cy[j], cz[k], out[idx]);
			}

	fclose(ptr);
	return 1;
}

template< typename T >
int nse::write_tecplot(const std::string& filename,
	const T* uout, const T* vout, const T* wout,
	const T* cx, const T* cy, const T* cz, // [array, coordinates]
	const int nx, const int ny, const int nz,
	const int ib, const int ie,
	const int jb, const int je,
	const int kb, const int ke,

	const char* title, const char* uname, const char* vname, const char* wname,
	const char* name_dimx, const char* name_dimy, const char* name_dimz,
	const T time)
{
	FILE* ptr = fopen(filename.c_str(), "w");
	if (ptr == NULL) return 0;

	fprintf(ptr, " TITLE = \"%s [F(%s,%s,%s)]\"\n",
		title, name_dimx, name_dimy, name_dimz);
	fprintf(ptr, " VARIABLES = \"%s\", \"%s\", \"%s\", \"%s\", \"%s\", \"%s\"\n",
		name_dimx, name_dimy, name_dimz, uname, vname, wname);
	fprintf(ptr, " ZONE I = %i, J = %i, K = %i, DATAPACKING = POINT, SOLUTIONTIME = %f\n",
		ie - ib + 1, je - jb + 1, ke - kb + 1, time);

	const int nyz = ny * nz;
	int i, j, k, idx;

	for (k = kb; k <= ke; k++)
		for (j = jb; j <= je; j++)
			for (i = ib; i <= ie; i++)
			{
				idx = i * nyz + j * nz + k;
				fprintf(ptr, "%f %f %f %f %f %f\n", cx[i], cy[j], cz[k],
					(T) 0.5 * (uout[idx] + uout[idx + nyz]),
					(T) 0.5 * (vout[idx] + vout[idx + nz]),
					(T) 0.5 * (wout[idx] + wout[idx + 1]));
			}

	fclose(ptr);
	return 1;
}
// -------------------------------------------------------------------------------------------- //

// Binary //
template< typename T >
int nse::write_binary_stamp(const std::string& filename,
	const binStamp< int >& index_stamp,
	const binStamp< double >& cpu_stamp,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId<T>& id, const T time)
{
	FILE* ptr = fopen(filename.c_str(), "wb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fwrite(id.header, sizeof(int), GridId<T>::hsize, ptr);	// header data //
	nstatus += fwrite(id.domain, sizeof(T), GridId<T>::dsize, ptr);	// domain definition //
	nstatus += fwrite(id.grid, sizeof(int), GridId<T>::gsize, ptr);	// grid definition //

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(cz, sizeof(T), nz, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, ptr);
	nstatus += fwrite(ez, sizeof(T), nz, ptr);

	// time stamp //
	T time_stamp = time;
	nstatus += fwrite(&time_stamp, sizeof(T), 1, ptr);

	nstatus += index_stamp.fwrite(ptr);	// index stamp //
	nstatus += cpu_stamp.fwrite(ptr);	// cpu stamp //


	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + index_stamp.size + cpu_stamp.size + 3;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::write_binary_stamp(const std::string& filename,
	const binNamedStamp< int >& index_stamp,
	const binNamedStamp< double >& cpu_stamp,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId<T>& id, const T time)
{
	FILE* ptr = fopen(filename.c_str(), "wb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fwrite(id.header, sizeof(int), GridId<T>::hsize, ptr);	// header data //
	nstatus += fwrite(id.domain, sizeof(T), GridId<T>::dsize, ptr);	// domain definition //
	nstatus += fwrite(id.grid, sizeof(int), GridId<T>::gsize, ptr);	// grid definition //

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(cz, sizeof(T), nz, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, ptr);
	nstatus += fwrite(ez, sizeof(T), nz, ptr);

	// time stamp //
	T time_stamp = time;
	nstatus += fwrite(&time_stamp, sizeof(T), 1, ptr);

	nstatus += index_stamp.fwrite(ptr);	// index stamp //
	nstatus += cpu_stamp.fwrite(ptr);	// cpu stamp //

	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + index_stamp.get_record_size() + cpu_stamp.get_record_size() + 1;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::write_binary(const std::string& filename,
	const T* xin, const char* name,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId<T>& id, const T time)
{
	FILE* ptr = fopen(filename.c_str(), "wb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fwrite(id.header, sizeof(int), GridId<T>::hsize, ptr);	// header data //
	nstatus += fwrite(id.domain, sizeof(T), GridId<T>::dsize, ptr);	// domain definition //
	nstatus += fwrite(id.grid, sizeof(int), GridId<T>::gsize, ptr);	// grid definition //

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(cz, sizeof(T), nz, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, ptr);
	nstatus += fwrite(ez, sizeof(T), nz, ptr);

	// time stamp //
	T time_stamp = time;
	nstatus += fwrite(&time_stamp, sizeof(T), 1, ptr);

	// field definition //
	int type = 0;               // scalar field
	int name_length = strlen(name);

	nstatus += fwrite(&type, sizeof(int), 1, ptr);
	nstatus += fwrite(&name_length, sizeof(int), 1, ptr);
	nstatus += fwrite(name, sizeof(char), name_length, ptr);

	// main data //
	nstatus += fwrite(xin, sizeof(T), nx * ny * nz, ptr);

	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + nx * ny * nz + name_length + 3;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::write_binary(const std::string& filename,
	const T* uin, const T* vin, const T* win,
	const char* uname, const char* vname, const char* wname,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId<T>& id, const T time)
{
	FILE* ptr = fopen(filename.c_str(), "wb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fwrite(id.header, sizeof(int), GridId<T>::hsize, ptr);	// header data //
	nstatus += fwrite(id.domain, sizeof(T), GridId<T>::dsize, ptr);	// domain definition //
	nstatus += fwrite(id.grid, sizeof(int), GridId<T>::gsize, ptr);	// grid definition //

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(cz, sizeof(T), nz, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, ptr);
	nstatus += fwrite(ez, sizeof(T), nz, ptr);

	// time stamp //
	T time_stamp = time;
	nstatus += fwrite(&time_stamp, sizeof(T), 1, ptr);

	// field definition //
	int type = 1;               // vector field
	int name_length[3];
	name_length[0] = strlen(uname);
	name_length[1] = strlen(vname);
	name_length[2] = strlen(wname);

	nstatus += fwrite(&type, sizeof(int), 1, ptr);
	nstatus += fwrite(name_length, sizeof(int), 3, ptr);
	nstatus += fwrite(uname, sizeof(char), name_length[0], ptr);
	nstatus += fwrite(vname, sizeof(char), name_length[1], ptr);
	nstatus += fwrite(wname, sizeof(char), name_length[2], ptr);

	// main data //
	nstatus += fwrite(uin, sizeof(T), nx * ny * nz, ptr);
	nstatus += fwrite(vin, sizeof(T), nx * ny * nz, ptr);
	nstatus += fwrite(win, sizeof(T), nx * ny * nz, ptr);

	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + 3 * nx * ny * nz +
		name_length[0] + name_length[1] + name_length[2] + 5;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::read_binary_stamp(const std::string& filename,
	binStamp< int >& index_stamp,
	binStamp< double >& cpu_stamp,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId<T>& id, T* time)
{
	FILE* ptr = fopen(filename.c_str(), "rb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fread(id.header, sizeof(int), GridId<T>::hsize_r3d, ptr);
	if (!id.check(3)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

	nstatus += fread_sp(id.domain, id.data_type_size(), GridId<T>::dsize_r3d, ptr);
	nstatus += fread(id.grid, sizeof(int), GridId<T>::gsize_r3d, ptr);

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	allocate(cx, cy, cz, nx, ny, nz);
	allocate(ex, ey, ez, nx, ny, nz);

#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*cz), id.data_type_size(), nz, ptr);
#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*ez), id.data_type_size(), nz, ptr);

	nstatus += fread_sp(time, id.data_type_size(), 1, ptr);	// time stamp

															// index stamp //
#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += index_stamp.fread(ptr);
#else
	nstatus += index_stamp.fread(ptr, 1);
#endif

	// cpu stamp //
	nstatus += cpu_stamp.fread(ptr);

	fclose(ptr);
	id.reset_data_type_size();	// re-setting data type 

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		index_stamp.size + cpu_stamp.size + 2 +
#ifndef _USE_DEPRECATED_WST_FORMAT
		2 * (nx + ny + nz) + 1;
#else
		2 * nz;
#endif
	if (nstatus == nstatus_check) return 1;

	deallocate((*cx), (*cy), (*cz));
	deallocate((*ex), (*ey), (*ez));
	return 0;
}

template< typename T >
int nse::read_binary_stamp(const std::string& filename,
	binNamedStamp< int >& index_stamp,
	binNamedStamp< double >& cpu_stamp,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId<T>& id, T* time)
{
	FILE* ptr = fopen(filename.c_str(), "rb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fread(id.header, sizeof(int), GridId<T>::hsize_r3d, ptr);
	if (!id.check(3)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

	nstatus += fread_sp(id.domain, id.data_type_size(), GridId<T>::dsize_r3d, ptr);
	nstatus += fread(id.grid, sizeof(int), GridId<T>::gsize_r3d, ptr);

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	allocate(cx, cy, cz, nx, ny, nz);
	allocate(ex, ey, ez, nx, ny, nz);

	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
	nstatus += fread_sp((*cz), id.data_type_size(), nz, ptr);

	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, ptr);
	nstatus += fread_sp((*ez), id.data_type_size(), nz, ptr);

	nstatus += fread_sp(time, id.data_type_size(), 1, ptr);	// time stamp

	// index stamp //
	nstatus += index_stamp.fread(ptr);
	// cpu stamp //
	nstatus += cpu_stamp.fread(ptr);

	fclose(ptr);
	id.reset_data_type_size();	// re-setting data type 

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		index_stamp.get_record_size() + cpu_stamp.get_record_size() +
		2 * (nx + ny + nz) + 1;
	if (nstatus == nstatus_check) return 1;

	deallocate((*cx), (*cy), (*cz));
	deallocate((*ex), (*ey), (*ez));
	return 0;
}

template< typename T >
int nse::read_binary(const std::string& filename,
	T** xout, char** name,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId<T>& id, T* time)
{
	FILE* ptr = fopen(filename.c_str(), "rb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fread(id.header, sizeof(int), GridId<T>::hsize_r3d, ptr);
	if (!id.check(3)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

	nstatus += fread_sp(id.domain, id.data_type_size(), GridId<T>::dsize_r3d, ptr);
	nstatus += fread(id.grid, sizeof(int), GridId<T>::gsize_r3d, ptr);

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	allocate(cx, cy, cz, nx, ny, nz);
	allocate(ex, ey, ez, nx, ny, nz);

#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*cz), id.data_type_size(), nz, ptr);
#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*ez), id.data_type_size(), nz, ptr);

	nstatus += fread_sp(time, id.data_type_size(), 1, ptr);	// time stamp

															// field definition //
	int field_type;
	nstatus += fread(&field_type, sizeof(int), 1, ptr);
	if (field_type != 0) {	// not scalar field //
		deallocate((*cx), (*cy), (*cz));
		deallocate((*ex), (*ey), (*ez));
		fclose(ptr);
		return 0;
	}

	int name_length;
	nstatus += fread(&name_length, sizeof(int), 1, ptr);

	(*name) = new char[name_length + 1];
	nstatus += fread((*name), sizeof(char), name_length, ptr);
	(*name)[name_length] = '\0';

	// main data //
	allocate(xout, nx * ny * nz);
	nstatus += fread_sp((*xout), id.data_type_size(), nx * ny * nz, ptr);

	fclose(ptr);
	id.reset_data_type_size();	// re-setting data type

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		nx * ny * nz + name_length + 3 +
#ifndef _USE_DEPRECATED_WST_FORMAT
		2 * (nx + ny + nz);
#else
		2 * nz;
#endif
	if (nstatus == nstatus_check) return 1;

	deallocate((*cx), (*cy), (*cz));
	deallocate((*ex), (*ey), (*ez));
	delete[](*name);
	deallocate((*xout));
	return 0;
}

template< typename T >
int nse::read_binary(const std::string& filename,
	T** uout, T** vout, T** wout,
	char** uname, char** vname, char** wname,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId<T>& id, T* time)
{
	FILE* ptr = fopen(filename.c_str(), "rb");
	if (ptr == NULL) return 0;

	int nstatus = 0;
	// header, domain & grid id //
	nstatus += fread(id.header, sizeof(int), GridId<T>::hsize_r3d, ptr);
	if (!id.check(3)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

	nstatus += fread_sp(id.domain, id.data_type_size(), GridId<T>::dsize_r3d, ptr);
	nstatus += fread(id.grid, sizeof(int), GridId<T>::gsize_r3d, ptr);

	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	// grid coordinates //
	allocate(cx, cy, cz, nx, ny, nz);
	allocate(ex, ey, ez, nx, ny, nz);

#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*cz), id.data_type_size(), nz, ptr);
#ifndef _USE_DEPRECATED_WST_FORMAT
	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, ptr);
#endif
	nstatus += fread_sp((*ez), id.data_type_size(), nz, ptr);

	nstatus += fread_sp(time, id.data_type_size(), 1, ptr);	// time stamp

															// field definition //
	int field_type;
	nstatus += fread(&field_type, sizeof(int), 1, ptr);
	if (field_type != 1) {	// not vector field //
		deallocate((*cx), (*cy), (*cz));
		deallocate((*ex), (*ey), (*ez));
		fclose(ptr);
		return 0;
	}

	int name_length[3];
	nstatus += fread(name_length, sizeof(int), 3, ptr);

	(*uname) = new char[name_length[0] + 1];
	(*vname) = new char[name_length[1] + 1];
	(*wname) = new char[name_length[2] + 1];
	nstatus += fread((*uname), sizeof(char), name_length[0], ptr);
	nstatus += fread((*vname), sizeof(char), name_length[1], ptr);
	nstatus += fread((*wname), sizeof(char), name_length[2], ptr);
	(*uname)[name_length[0]] = '\0';
	(*vname)[name_length[1]] = '\0';
	(*wname)[name_length[2]] = '\0';

	// main data //
	allocate(uout, vout, wout, nx * ny * nz);
	nstatus += fread_sp((*uout), id.data_type_size(), nx * ny * nz, ptr);
	nstatus += fread_sp((*vout), id.data_type_size(), nx * ny * nz, ptr);
	nstatus += fread_sp((*wout), id.data_type_size(), nx * ny * nz, ptr);

	fclose(ptr);
	id.reset_data_type_size();	// re-setting data type

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		3 * nx * ny * nz + name_length[0] + name_length[1] + name_length[2] + 5 +
#ifndef _USE_DEPRECATED_WST_FORMAT
		2 * (nx + ny + nz);
#else
		2 * nz;
#endif
	if (nstatus == nstatus_check) return 1;

	deallocate((*cx), (*cy), (*cz));
	deallocate((*ex), (*ey), (*ez));
	delete[](*uname); delete[](*vname); delete[](*wname);
	deallocate((*uout), (*vout), (*wout));
	return 0;
}
// -------------------------------------------------------------------------------------------- //

// MPI-I/O 3D datatype //
template< typename T >
void nse::mpi_io_write_datatype(MPI_Datatype* file_view, MPI_Datatype* local_view,

	const int mpi_nx, const int mpi_ny, const int mpi_nz,
	const int nx, const int ny, const int nz,
	const int gcx, const int gcy, const int gcz,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z)
{
	const int mpi_dim_size[3] = { mpi_nx, mpi_ny, mpi_nz };
	const int dim_size[3] = { nx, ny, nz };
	const int nghost[3] = { gcx, gcy, gcz };
	const MPI_Comm comm[3] = { comm_x, comm_y, comm_z };

	mpi_io_write_datatype< T, 3 >(file_view, local_view,
		mpi_dim_size, dim_size, nghost, comm);
}

template< typename T >
void nse::mpi_io_read_datatype(MPI_Datatype* file_view, MPI_Datatype* local_view,

	const int mpi_nx, const int mpi_ny, const int mpi_nz,
	const int nx, const int ny, const int nz,
	const int gcx, const int gcy, const int gcz,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z)
{
	const int mpi_dim_size[3] = { mpi_nx, mpi_ny, mpi_nz };
	const int dim_size[3] = { nx, ny, nz };
	const int nghost[3] = { gcx, gcy, gcz };
	const MPI_Comm comm[3] = { comm_x, comm_y, comm_z };

	mpi_io_read_datatype< T, 3 >(file_view, local_view,
		mpi_dim_size, dim_size, nghost, comm);
}
// -------------------------------------------------------------------------------------------- //

// MPI-Binary //
template< typename T >
int nse::mpi_write_binary(const std::string& filename,
	const char* mpi_datarep,
	const MPI_Comm comm, const int header_rank,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

	const T* xin, const char* name,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId< T >& id, const T time)
{
	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	MPI_Offset header_size = GridId< T >::id_byte_size +
		2 * (nx + ny + nz) * sizeof(T) +
		strlen(name) * sizeof(char) +
		2 * sizeof(int) + sizeof(T);

	MPI_File ptr;
	int status = MPI_File_open(comm, (char*)filename.c_str(),
		MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &ptr);
	if (status != MPI_SUCCESS) return 0;	// MPI file open failure

	int nstatus = 0;
	int rank;
	MPI_Comm_rank(comm, &rank);
	if (rank == header_rank) {	// header

								// header, domain & grid id //
		status = MPI_File_write(ptr, (void*)id.header, GridId< T >::hsize, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::hsize;
		status = MPI_File_write(ptr, (void*)id.domain, GridId< T >::dsize, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::dsize;
		status = MPI_File_write(ptr, (void*)id.grid, GridId< T >::gsize, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::gsize;

		// grid coordinates //
		status = MPI_File_write(ptr, (void*)cx, nx, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nx;
		status = MPI_File_write(ptr, (void*)cy, ny, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += ny;
		status = MPI_File_write(ptr, (void*)cz, nz, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nz;
		status = MPI_File_write(ptr, (void*)ex, nx, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nx;
		status = MPI_File_write(ptr, (void*)ey, ny, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += ny;
		status = MPI_File_write(ptr, (void*)ez, nz, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nz;

		// time stamp //
		T time_stamp = time;
		status = MPI_File_write(ptr, &time_stamp, 1, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 1;

		// field definition //
		int type = 0;		// scalar field
		int name_length = strlen(name);

		status = MPI_File_write(ptr, &type, 1, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 1;
		status = MPI_File_write(ptr, &name_length, 1, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 1;
		status = MPI_File_write(ptr, (void*)name, name_length, MPI_CHAR, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += name_length;
	}
	MPI_File_sync(ptr);
	mpi_broadcast(&nstatus, 1, header_rank, comm);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);
	int pnz = par_local_size_comm(nz, gcz, comm_z);

	// main data description //
	MPI_Datatype file_view, local_view;
	mpi_io_write_datatype< T >(&file_view, &local_view,
		nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
		comm_x, comm_y, comm_z);

	// main data //
	MPI_File_set_view(ptr, header_size, mpi_type< T >(),
		file_view, (char*)mpi_datarep, MPI_INFO_NULL);

	status = MPI_File_write_all(ptr, (void*)xin, 1, local_view, MPI_STATUS_IGNORE);
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;

	MPI_File_close(&ptr);
	MPI_Type_free(&file_view);
	MPI_Type_free(&local_view);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + nx * ny * nz + strlen(name) + 3;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::mpi_write_binary(const std::string& filename,
	const char* mpi_datarep,
	const MPI_Comm comm, const int header_rank,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

	const T* uin, const T* vin, const T* win,
	const char* uname, const char* vname, const char* wname,

	const T* cx, const T* cy, const T* cz,
	const T* ex, const T* ey, const T* ez,
	const GridId< T >& id, const T time)
{
	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	MPI_Offset header_size = GridId< T >::id_byte_size +
		2 * (nx + ny + nz) * sizeof(T) +
		(strlen(uname) + strlen(vname) + strlen(wname)) * sizeof(char) +
		4 * sizeof(int) + sizeof(T);

	MPI_File ptr;
	int status = MPI_File_open(comm, (char*)filename.c_str(),
		MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &ptr);
	if (status != MPI_SUCCESS) return 0;	// MPI file open failure

	int nstatus = 0;
	int rank;
	MPI_Comm_rank(comm, &rank);
	if (rank == header_rank) {	// header

								// header, domain & grid id //
		status = MPI_File_write(ptr, (void*)id.header, GridId< T >::hsize, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::hsize;
		status = MPI_File_write(ptr, (void*)id.domain, GridId< T >::dsize, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::dsize;
		status = MPI_File_write(ptr, (void*)id.grid, GridId< T >::gsize, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId< T >::gsize;

		// grid coordinates //
		status = MPI_File_write(ptr, (void*)cx, nx, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nx;
		status = MPI_File_write(ptr, (void*)cy, ny, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += ny;
		status = MPI_File_write(ptr, (void*)cz, nz, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nz;
		status = MPI_File_write(ptr, (void*)ex, nx, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nx;
		status = MPI_File_write(ptr, (void*)ey, ny, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += ny;
		status = MPI_File_write(ptr, (void*)ez, nz, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += nz;

		// time stamp //
		T time_stamp = time;
		status = MPI_File_write(ptr, &time_stamp, 1, mpi_type< T >(), MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 1;

		// field definition //
		int type = 1;		// vector field
		int name_length[3];
		name_length[0] = strlen(uname);
		name_length[1] = strlen(vname);
		name_length[2] = strlen(wname);

		status = MPI_File_write(ptr, &type, 1, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 1;
		status = MPI_File_write(ptr, name_length, 3, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 3;
		status = MPI_File_write(ptr, (void*)uname, name_length[0], MPI_CHAR, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += name_length[0];
		status = MPI_File_write(ptr, (void*)vname, name_length[1], MPI_CHAR, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += name_length[1];
		status = MPI_File_write(ptr, (void*)wname, name_length[2], MPI_CHAR, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += name_length[2];
	}
	MPI_File_sync(ptr);
	mpi_broadcast(&nstatus, 1, header_rank, comm);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);
	int pnz = par_local_size_comm(nz, gcz, comm_z);

	// main data description //
	MPI_Datatype file_view, local_view;
	mpi_io_write_datatype< T >(&file_view, &local_view,
		nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
		comm_x, comm_y, comm_z);

	// main data //
	MPI_File_set_view(ptr, header_size, mpi_type< T >(),
		file_view, (char*)mpi_datarep, MPI_INFO_NULL);

	status = MPI_File_write_all(ptr, (void*)uin, 1, local_view, MPI_STATUS_IGNORE);
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;
	status = MPI_File_write_all(ptr, (void*)vin, 1, local_view, MPI_STATUS_IGNORE);
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;
	status = MPI_File_write_all(ptr, (void*)win, 1, local_view, MPI_STATUS_IGNORE);
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;

	MPI_File_close(&ptr);
	MPI_Type_free(&file_view);
	MPI_Type_free(&local_view);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny + nz) + 3 * nx * ny * nz +
		strlen(uname) + strlen(vname) + strlen(wname) + 5;
	return (nstatus == nstatus_check);
}

template< typename T >
int nse::mpi_read_binary(const std::string& filename,
	const char* mpi_datarep,
	const MPI_Comm comm, const int header_rank,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

	T** xout, char** name,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId< T >& id, T* time)
{
	MPI_Offset header_size = GridId<T>::id_byte_size_r3d +
		2 * sizeof(int);

	MPI_File ptr;
	int status = MPI_File_open(comm, (char*)filename.c_str(),
		MPI_MODE_RDONLY, MPI_INFO_NULL, &ptr);
	if (status != MPI_SUCCESS) return 0;	// MPI file open failure

	int name_length;
	int nstatus = 0;
	int rank, header_offset = 0, status_id = 0;
	MPI_Comm_rank(comm, &rank);
	if (rank == header_rank) {

		// header, domain & grid id //
		status = MPI_File_read(ptr, id.header, GridId<T>::hsize_r3d, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId<T>::hsize_r3d;
		if (id.check(3)) {	// check id //
			status = mpi_fread_sp(ptr, id.domain, GridId<T>::dsize_r3d, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += GridId<T>::dsize_r3d;
			status = MPI_File_read(ptr, id.grid, GridId<T>::gsize_r3d, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += GridId<T>::gsize_r3d;

			header_offset += GridId<T>::dsize_r3d *
				(id.data_type_size() - sizeof(T));	// correcting header size due to grid id

													// grid parameters //
			int nx, ny, nz, gcx, gcy, gcz;
			id.grid_dim(1, &nx, &gcx);
			id.grid_dim(2, &ny, &gcy);
			id.grid_dim(3, &nz, &gcz);

			// grid coordinates //
			allocate(cx, cy, cz, nx, ny, nz);
			allocate(ex, ey, ez, nx, ny, nz);

#ifndef _USE_DEPRECATED_WST_FORMAT
			status = mpi_fread_sp(ptr, (*cx), nx, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nx;
			status = mpi_fread_sp(ptr, (*cy), ny, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += ny;
#endif
			status = mpi_fread_sp(ptr, (*cz), nz, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nz;
#ifndef _USE_DEPRECATED_WST_FORMAT
			status = mpi_fread_sp(ptr, (*ex), nx, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nx;
			status = mpi_fread_sp(ptr, (*ey), ny, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += ny;
#endif
			status = mpi_fread_sp(ptr, (*ez), nz, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nz;
#ifndef _USE_DEPRECATED_WST_FORMAT
			header_offset += 2 * (nx + ny + nz) * id.data_type_size();
#else
			header_offset += 2 * nz * id.data_type_size();
#endif

			// time stamp //
			status = mpi_fread_sp(ptr, time, 1, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += 1;
			header_offset += id.data_type_size();

			// field definition //
			int field_type;
			status = MPI_File_read(ptr, &field_type, 1, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += 1;
			if (field_type == 0)	// scalar field //
			{
				status = MPI_File_read(ptr, &name_length, 1, MPI_INT, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += 1;
				header_offset += name_length * sizeof(char);

				(*name) = new char[name_length + 1];
				status = MPI_File_read(ptr, (*name), name_length, MPI_CHAR, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += name_length;
				(*name)[name_length] = '\0';

				status_id = 1;
			}
			else
			{
				deallocate((*cx), (*cy), (*cz));
				deallocate((*ex), (*ey), (*ez));
			}
		}
	}
	mpi_broadcast(&status_id, 1, header_rank, comm);

	if (!status_id) {
		MPI_File_close(&ptr);
		return 0;
	}
	// read status - OK - //
	id.mpi_broadcast(header_rank, comm);

	mpi_broadcast(&nstatus, 1, header_rank, comm);
	mpi_broadcast(&name_length, 1, header_rank, comm);

	// correct header size //
	mpi_broadcast(&header_offset, 1, header_rank, comm);
	header_size += (MPI_Offset)header_offset;

	// main data description //
	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);
	int pnz = par_local_size_comm(nz, gcz, comm_z);

	MPI_Datatype file_view, local_view;
	if (id.data_type_size() == sizeof(float)) {	// input=[float]
		mpi_io_read_datatype< float >(&file_view, &local_view,
			nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
			comm_x, comm_y, comm_z);

		MPI_File_set_view(ptr, header_size, mpi_type< float >(),
			file_view, (char*)mpi_datarep, MPI_INFO_NULL);
	}
	if (id.data_type_size() == sizeof(double)) {	// input=[double]
		mpi_io_read_datatype< double >(&file_view, &local_view,
			nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
			comm_x, comm_y, comm_z);

		MPI_File_set_view(ptr, header_size, mpi_type< double >(),
			file_view, (char*)mpi_datarep, MPI_INFO_NULL);
	}

	// main data //
	allocate(xout, pnx * pny * pnz);
	status = mpi_fread_all_sp(ptr, (*xout), pnx * pny * pnz,
		local_view, id.data_type_size());
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;

	MPI_File_close(&ptr);
	MPI_Type_free(&file_view);
	MPI_Type_free(&local_view);

	id.reset_data_type_size();	// re-setting data type

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		nx * ny * nz + name_length + 3 +
#ifndef _USE_DEPRECATED_WST_FORMAT
		2 * (nx + ny + nz);
#else
		2 * nz;
#endif
	if (nstatus == nstatus_check) return 1;

	if (rank == header_rank) {
		deallocate((*cx), (*cy), (*cz));
		deallocate((*ex), (*ey), (*ez));
		delete[](*name);
	}
	deallocate((*xout));
	return 0;
}

template< typename T >
int nse::mpi_read_binary(const std::string& filename,
	const char* mpi_datarep,
	const MPI_Comm comm, const int header_rank,
	const MPI_Comm comm_x, const MPI_Comm comm_y, const MPI_Comm comm_z,

	T** uout, T** vout, T** wout,
	char** uname, char** vname, char** wname,

	T** cx, T** cy, T** cz, T** ex, T** ey, T** ez,
	GridId< T >& id, T* time)
{
	MPI_Offset header_size = GridId<T>::id_byte_size_r3d +
		4 * sizeof(int);

	MPI_File ptr;
	int status = MPI_File_open(comm, (char*)filename.c_str(),
		MPI_MODE_RDONLY, MPI_INFO_NULL, &ptr);
	if (status != MPI_SUCCESS) return 0;	// MPI file open failure

	int name_length[3];
	int nstatus = 0;
	int rank, header_offset = 0, status_id = 0;
	MPI_Comm_rank(comm, &rank);
	if (rank == header_rank) {

		// header, domain & grid id //
		status = MPI_File_read(ptr, id.header, GridId<T>::hsize_r3d, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId<T>::hsize_r3d;
		if (id.check(3)) {	// check id //
			status = mpi_fread_sp(ptr, id.domain, GridId<T>::dsize_r3d, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += GridId<T>::dsize_r3d;
			status = MPI_File_read(ptr, id.grid, GridId<T>::gsize_r3d, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += GridId<T>::gsize_r3d;

			header_offset += GridId<T>::dsize_r3d *
				(id.data_type_size() - sizeof(T));	// correcting header size due to grid id

													// grid parameters //
			int nx, ny, nz, gcx, gcy, gcz;
			id.grid_dim(1, &nx, &gcx);
			id.grid_dim(2, &ny, &gcy);
			id.grid_dim(3, &nz, &gcz);

			// grid coordinates //
			allocate(cx, cy, cz, nx, ny, nz);
			allocate(ex, ey, ez, nx, ny, nz);

#ifndef _USE_DEPRECATED_WST_FORMAT
			status = mpi_fread_sp(ptr, (*cx), nx, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nx;
			status = mpi_fread_sp(ptr, (*cy), ny, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += ny;
#endif
			status = mpi_fread_sp(ptr, (*cz), nz, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nz;
#ifndef _USE_DEPRECATED_WST_FORMAT
			status = mpi_fread_sp(ptr, (*ex), nx, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nx;
			status = mpi_fread_sp(ptr, (*ey), ny, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += ny;
#endif
			status = mpi_fread_sp(ptr, (*ez), nz, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += nz;
#ifndef _USE_DEPRECATED_WST_FORMAT
			header_offset += 2 * (nx + ny + nz) * id.data_type_size();
#else
			header_offset += 2 * nz * id.data_type_size();
#endif

			// time stamp //
			status = mpi_fread_sp(ptr, time, 1, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += 1;
			header_offset += id.data_type_size();

			// field definition //
			int field_type;
			status = MPI_File_read(ptr, &field_type, 1, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += 1;
			if (field_type == 1)	// vector field //
			{
				status = MPI_File_read(ptr, name_length, 3, MPI_INT, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += 3;
				header_offset += sizeof(char) *
					(name_length[0] + name_length[1] + name_length[2]);

				(*uname) = new char[name_length[0] + 1];
				(*vname) = new char[name_length[1] + 1];
				(*wname) = new char[name_length[2] + 1];
				status = MPI_File_read(ptr, (*uname), name_length[0], MPI_CHAR, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += name_length[0];
				status = MPI_File_read(ptr, (*vname), name_length[1], MPI_CHAR, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += name_length[1];
				status = MPI_File_read(ptr, (*wname), name_length[2], MPI_CHAR, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += name_length[2];
				(*uname)[name_length[0]] = '\0';
				(*vname)[name_length[1]] = '\0';
				(*wname)[name_length[2]] = '\0';

				status_id = 1;
			}
			else
			{
				deallocate((*cx), (*cy), (*cz));
				deallocate((*ex), (*ey), (*ez));
			}
		}
	}
	mpi_broadcast(&status_id, 1, header_rank, comm);

	if (!status_id) {
		MPI_File_close(&ptr);
		return 0;
	}
	// read status - OK - //
	id.mpi_broadcast(header_rank, comm);

	mpi_broadcast(&nstatus, 1, header_rank, comm);
	mpi_broadcast(name_length, 3, header_rank, comm);

	// correct header size //
	mpi_broadcast(&header_offset, 1, header_rank, comm);
	header_size += (MPI_Offset)header_offset;

	// main data description //
	int nx, ny, nz, gcx, gcy, gcz;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);
	id.grid_dim(3, &nz, &gcz);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);
	int pnz = par_local_size_comm(nz, gcz, comm_z);

	MPI_Datatype file_view, local_view;
	if (id.data_type_size() == sizeof(float)) {	// input=[float]
		mpi_io_read_datatype< float >(&file_view, &local_view,
			nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
			comm_x, comm_y, comm_z);

		MPI_File_set_view(ptr, header_size, mpi_type< float >(),
			file_view, (char*)mpi_datarep, MPI_INFO_NULL);
	}
	if (id.data_type_size() == sizeof(double)) {	// input=[double]
		mpi_io_read_datatype< double >(&file_view, &local_view,
			nx, ny, nz, pnx, pny, pnz, gcx, gcy, gcz,
			comm_x, comm_y, comm_z);

		MPI_File_set_view(ptr, header_size, mpi_type< double >(),
			file_view, (char*)mpi_datarep, MPI_INFO_NULL);
	}

	// main data //
	allocate(uout, vout, wout, pnx * pny * pnz);
	status = mpi_fread_all_sp(ptr, (*uout), pnx * pny * pnz,
		local_view, id.data_type_size());
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;
	status = mpi_fread_all_sp(ptr, (*vout), pnx * pny * pnz,
		local_view, id.data_type_size());
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;
	status = mpi_fread_all_sp(ptr, (*wout), pnx * pny * pnz,
		local_view, id.data_type_size());
	if (status == MPI_SUCCESS) nstatus += nx * ny * nz;

	MPI_File_close(&ptr);
	MPI_Type_free(&file_view);
	MPI_Type_free(&local_view);

	id.reset_data_type_size();	// re-setting data type

	const int nstatus_check =
		GridId<T>::hsize_r3d + GridId<T>::dsize_r3d + GridId<T>::gsize_r3d +
		3 * nx * ny * nz + name_length[0] + name_length[1] + name_length[2] + 5 +
#ifndef _USE_DEPRECATED_WST_FORMAT
		2 * (nx + ny + nz);
#else
		2 * nz;
#endif
	if (nstatus == nstatus_check) return 1;

	if (rank == header_rank) {
		deallocate((*cx), (*cy), (*cz));
		deallocate((*ex), (*ey), (*ez));
		delete[](*uname); delete[](*vname); delete[](*wname);
	}
	deallocate((*uout), (*vout), (*wout));
	return 0;
}
// -------------------------------------------------------------------------------------------- //