#pragma once

// [io-base2d.h]: I/O: 2D 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"

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

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

	template< typename T >
	int write_tecplot(const std::string& filename,
		T** out, const T* cx, const T* cy, const int nvar,
		const int nx, const int ny,
		const int ib, const int ie,
		const int jb, const int je,

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

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

		const char* title, const char* uname, const char* vname,
		const char* name_dimx, const char* name_dimy,
		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* ex, const T* ey,
		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* ex, const T* ey,
		const GridId<T>& id, const T time);

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

		const T* cx, const T* cy,
		const T* ex, const T* ey,
		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** ex, T** ey,
		GridId<T>& id, T* time);

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

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

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

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

	// MPI-I/O 2D 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 nx, const int ny,
		const int gcx, const int gcy,
		const MPI_Comm comm_x, const MPI_Comm comm_y);

	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 nx, const int ny,
		const int gcx, const int gcy,
		const MPI_Comm comm_x, const MPI_Comm comm_y);
	// -------------------------------------------------------------------------------------------- //

	// 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 T* xin, const char* name,

		const T* cx, const T* cy,
		const T* ex, const T* ey,
		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 T* uin, const T* vin,
		const char* uname, const char* vname,

		const T* cx, const T* cy,
		const T* ex, const T* ey,
		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,

		T** xout, char** name,

		T** cx, T** cy, T** ex, T** ey,
		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,

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

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


	template< typename T >
	bool read_plain_2d(const std::string& filename,
		T** F, int* nx, int* ny);

	template< typename T >
	bool mpi_read_plain_2d(const std::string& filename,
		T** F, int* nx, int* ny, const MPI_Comm comm);
	// -------------------------------------------------------------------------------------------- //
}

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

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

	const char* title, const char* name,
	const char* name_dimx, const char* name_dimy,
	const T time)
{
	if ((ib < 0) || (jb < 0) || (ie >= nx) || (je >= ny)) return 0;

	FILE* ptr = fopen(filename.c_str(), "w");
	if (ptr == NULL) return 0;

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

	int i, j, idx;
	for (j = jb; j <= je; j++)
		for (i = ib; i <= ie; i++)
		{
			idx = i * ny + j;
			fprintf(ptr, "%f %f %f\n", cx[i], cy[j], out[idx]);
		}

	fclose(ptr);
	return 1;
}

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

	const char* title, const char** name,
	const char* name_dimx, const char* name_dimy,
	const T time)
{
	if ((ib < 0) || (jb < 0) || (ie >= nx) || (je >= ny)) return 0;

	FILE* ptr = fopen(filename.c_str(), "w");
	if (ptr == NULL) return 0;

	fprintf(ptr, " TITLE = \"%s [F_i(%s,%s), i=1 .. %i]\"\n", title, name_dimx, name_dimy, nvar);
	fprintf(ptr, " VARIABLES = \"%s\", \"%s\",", name_dimx, name_dimy);
	for (int k = 0; k < nvar - 1; k++)
		fprintf(ptr, " \"%s\",", name[k]);
	fprintf(ptr, " \"%s\"\n", name[nvar - 1]);
	fprintf(ptr, " ZONE I = %i, J = %i, DATAPACKING = POINT, SOLUTIONTIME = %f\n",
		ie - ib + 1, je - jb + 1, time);

	int i, j, idx;
	for (j = jb; j <= je; j++)
		for (i = ib; i <= ie; i++)
		{
			idx = i * ny + j;
			fprintf(ptr, "%f %f", cx[i], cy[j]);
			for (int k = 0; k < nvar; k++)
				fprintf(ptr, " %f", out[k][idx]);
			fprintf(ptr, "\n");
		}

	fclose(ptr);
	return 1;
}

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

	const char* title, const char* uname, const char* vname,
	const char* name_dimx, const char* name_dimy,
	const T time)
{
	if ((ib < 0) || (jb < 0) || (ie >= nx - 1) || (je >= ny - 1)) return 0;

	FILE* ptr = fopen(filename.c_str(), "w");
	if (ptr == NULL) return 0;

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

	int i, j, idx;
	for (j = jb; j <= je; j++)
		for (i = ib; i <= ie; i++)
		{
			idx = i * ny + j;
			fprintf(ptr, "%f %f %f %f\n", cx[i], cy[j],
				(T) 0.5 * (uout[idx] + uout[idx + ny]),
				(T) 0.5 * (vout[idx] + vout[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* ex, const T* ey,
	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, gcx, gcy;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, 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) + index_stamp.size + cpu_stamp.size + 3;
	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* ex, const T* ey,
	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, gcx, gcy;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, 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, ptr);

	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny) + nx * ny + 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 char* uname, const char* vname,

	const T* cx, const T* cy,
	const T* ex, const T* ey,
	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, gcx, gcy;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);

	// grid coordinates //
	nstatus += fwrite(cx, sizeof(T), nx, ptr);
	nstatus += fwrite(cy, sizeof(T), ny, ptr);
	nstatus += fwrite(ex, sizeof(T), nx, ptr);
	nstatus += fwrite(ey, sizeof(T), ny, 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[2];
	name_length[0] = strlen(uname);
	name_length[1] = strlen(vname);

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

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

	fclose(ptr);

	const int nstatus_check = GridId<T>::hsize + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny) + 2 * nx * ny +
		name_length[0] + name_length[1] + 4;
	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** ex, T** ey,
	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, ptr);
	if (!id.check(2)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

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

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

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

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

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

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

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

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

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

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

	T** cx, T** cy, T** ex, T** ey,
	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, ptr);
	if (!id.check(2)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

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

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

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

	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, 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));
		deallocate((*ex), (*ey));
		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);
	nstatus += fread_sp((*xout), id.data_type_size(), nx * ny, ptr);

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

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

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

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

	T** cx, T** cy, T** ex, T** ey,
	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, ptr);
	if (!id.check(2)) {	// check id failed //
		fclose(ptr);
		return 0;
	}

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

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

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

	nstatus += fread_sp((*cx), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*cy), id.data_type_size(), ny, ptr);
	nstatus += fread_sp((*ex), id.data_type_size(), nx, ptr);
	nstatus += fread_sp((*ey), id.data_type_size(), ny, 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));
		deallocate((*ex), (*ey));
		fclose(ptr);
		return 0;
	}

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

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

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

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

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

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

// MPI-I/O 2D 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 nx, const int ny,
	const int gcx, const int gcy,
	const MPI_Comm comm_x, const MPI_Comm comm_y)
{
	const int mpi_dim_size[2] = { mpi_nx, mpi_ny };
	const int dim_size[2] = { nx, ny };
	const int nghost[2] = { gcx, gcy };
	const MPI_Comm comm[2] = { comm_x, comm_y };

	mpi_io_write_datatype< T, 2 >(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 nx, const int ny,
	const int gcx, const int gcy,
	const MPI_Comm comm_x, const MPI_Comm comm_y)
{
	const int mpi_dim_size[2] = { mpi_nx, mpi_ny };
	const int dim_size[2] = { nx, ny };
	const int nghost[2] = { gcx, gcy };
	const MPI_Comm comm[2] = { comm_x, comm_y };

	mpi_io_read_datatype< T, 2 >(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 T* xin, const char* name,

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

	MPI_Offset header_size = GridId< T >::id_byte_size +
		2 * (nx + ny) * 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*)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;

		// 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);

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

	// 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;

	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) + nx * ny + 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 T* uin, const T* vin,
	const char* uname, const char* vname,

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

	MPI_Offset header_size = GridId< T >::id_byte_size +
		2 * (nx + ny) * sizeof(T) +
		(strlen(uname) + strlen(vname)) * sizeof(char) +
		3 * 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*)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;

		// 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[2];
		name_length[0] = strlen(uname);
		name_length[1] = strlen(vname);

		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, 2, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += 2;
		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];
	}
	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);

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

	// 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;
	status = MPI_File_write_all(ptr, (void*)vin, 1, local_view, MPI_STATUS_IGNORE);
	if (status == MPI_SUCCESS) nstatus += nx * ny;

	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) + 2 * nx * ny +
		strlen(uname) + strlen(vname) + 4;
	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,

	T** xout, char** name,

	T** cx, T** cy, T** ex, T** ey,
	GridId< T >& id, T* time)
{
	MPI_Offset header_size = GridId<T>::id_byte_size +
		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, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId<T>::hsize;
		if (id.check(2)) {	// check id //
			status = mpi_fread_sp(ptr, id.domain, GridId<T>::dsize, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += GridId<T>::dsize;
			status = MPI_File_read(ptr, id.grid, GridId<T>::gsize, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += GridId<T>::gsize;

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

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

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

			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;
			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;
			header_offset += 2 * (nx + ny) * id.data_type_size();

			// 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));
				deallocate((*ex), (*ey));
			}
		}
	}
	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, gcx, gcy;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);

	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, pnx, pny, gcx, gcy,
			comm_x, comm_y);

		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, pnx, pny, gcx, gcy,
			comm_x, comm_y);

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

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

	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 + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny) + nx * ny + name_length + 3;
	if (nstatus == nstatus_check) return 1;

	if (rank == header_rank) {
		deallocate((*cx), (*cy));
		deallocate((*ex), (*ey));
		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,

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

	T** cx, T** cy, T** ex, T** ey,
	GridId< T >& id, T* time)
{
	MPI_Offset header_size = GridId<T>::id_byte_size +
		3 * 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[2];
	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, MPI_INT, MPI_STATUS_IGNORE);
		if (status == MPI_SUCCESS) nstatus += GridId<T>::hsize;
		if (id.check(2)) {	// check id //
			status = mpi_fread_sp(ptr, id.domain, GridId<T>::dsize, id.data_type_size());
			if (status == MPI_SUCCESS) nstatus += GridId<T>::dsize;
			status = MPI_File_read(ptr, id.grid, GridId<T>::gsize, MPI_INT, MPI_STATUS_IGNORE);
			if (status == MPI_SUCCESS) nstatus += GridId<T>::gsize;

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

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

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

			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;
			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;
			header_offset += 2 * (nx + ny) * id.data_type_size();

			// 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, 2, MPI_INT, MPI_STATUS_IGNORE);
				if (status == MPI_SUCCESS) nstatus += 2;
				header_offset += sizeof(char)*
					(name_length[0] + name_length[1]);

				(*uname) = new char[name_length[0] + 1];
				(*vname) = new char[name_length[1] + 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];
				(*uname)[name_length[0]] = '\0';
				(*vname)[name_length[1]] = '\0';

				status_id = 1;
			}
			else
			{
				deallocate((*cx), (*cy));
				deallocate((*ex), (*ey));
			}
		}
	}
	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, 2, 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, gcx, gcy;
	id.grid_dim(1, &nx, &gcx);
	id.grid_dim(2, &ny, &gcy);

	int pnx = par_local_size_comm(nx, gcx, comm_x);
	int pny = par_local_size_comm(ny, gcy, comm_y);

	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, pnx, pny, gcx, gcy,
			comm_x, comm_y);

		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, pnx, pny, gcx, gcy,
			comm_x, comm_y);

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

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

	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 + GridId<T>::dsize + GridId<T>::gsize +
		2 * (nx + ny) + 2 * nx * ny +
		name_length[0] + name_length[1] + 4;
	if (nstatus == nstatus_check) return 1;

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

template< typename T >
bool nse::read_plain_2d(const std::string& filename,
	T** F, int* nx, int* ny)
{
	FILE *ptr = fopen(filename.c_str(), "rt");
	if (ptr == NULL) return false;

	if (fscanf(ptr, c_io_fmt< int >(), nx) != 1) {
		fclose(ptr);
		return false;
	}
	if (fscanf(ptr, c_io_fmt< int >(), ny) != 1) {
		fclose(ptr);
		return false;
	}

	if (((*nx) <= 0) || ((*ny) <= 0)) {
		fclose(ptr);
		return false;
	}

	allocate(F, (*nx) * (*ny));

	for (int k = 0; k < (*nx) * (*ny); k++)
	{
		if (fscanf(ptr, c_io_fmt< T >(), &((*F)[k])) != 1) {
			fclose(ptr);
			deallocate((*F));
			return false;
		}
	}

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

template< typename T >
bool nse::mpi_read_plain_2d(const std::string& filename,
	T** F, int* nx, int* ny, const MPI_Comm comm)
{
	bool status;

	int rank;
	MPI_Comm_rank(comm, &rank);
	if (rank == 0) {
		status = read_plain_2d(filename, F, nx, ny);
	}

	mpi_broadcast(&status, 1, 0, comm);
	if (!status) return false;

	mpi_broadcast(nx, 1, 0, comm);
	mpi_broadcast(ny, 1, 0, comm);
	if (rank != 0)
		allocate_vnull(F, (*nx) * (*ny));

	mpi_broadcast(*F, (*nx) * (*ny), 0, comm);
	return true;
}
// -------------------------------------------------------------------------------------------- //