/******************************************************************************
 * @internal
 * @file        rCFD_parallel.h
 * @copyright   (C) 2021
 *              Particulate Flow Modelling,
 *              Johannes Kepler University, Linz, Austria
 * @website     www.particulate-flow.at
 * @authors     Stefan Pirker
 *****************************************************************************/

#ifndef RCFD_PARALLEL
#define RCFD_PARALLEL

#include "rCFD.h"
#include "rCFD_types.h"
#include "rCFD_macros.h"

#if 1 /* macros */

    #define     _MPI_Cell               rcfd->Parallel->MPI_Topo->MPI_Cell[i_layer]

#endif

#if 1 /* typedef */
#if RP_NODE

    typedef struct MPI_cells_struct
    {
        int     total_number_of_ext_cells;          // node-0

        int     number_of_host_cells;

        int     *number_of_ext_cells_per_node;

        int     *number_of_host_cells_per_node;

        int     *host_of_cell;                      /* all nodes: list of host of all cells */

        int     *host_of_ext_cells;                 /* node-0: total list of all ext. cells */

        int     *hosting_cell_index;                /* all nodes: list of cells, which host ext. cells of other nodes */

        int     *host2ext_index;

        double  *data;                              // all nodes: cell_data of corona cells

    } MPI_cells_type;

    typedef struct MPI_topo_struct
    {
        MPI_cells_type  *MPI_Cell;

    } MPI_topo_type;

    typedef struct C2C_MPI_shift_struct
    {
        int         c0,node0;
        int         c1,node1;
        double      w0;

        double      *data;

    } C2C_MPI_shift_type;

    typedef struct C2C_MPI_struct
    {
        short               format;

        int                 max_number_of_MPI_shifts;

        C2C_MPI_shift_type  *shifts_in;

        int                 number_of_shifts_in; /* temporally use in C2Cs_prepare */

     } C2C_MPI_type;

#endif
#endif


#define TAG_GATHERV 101
#if defined(RCFD_STANDALONE)
#if defined(RCFD_LIB_MPI)
    int     prf_gisum1  (void *handle, int a)     { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT,       MPI_SUM, rcfd->Parallel->world); return a; }
    int     prf_gihigh1 (void *handle, int a)     { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT,       MPI_MAX, rcfd->Parallel->world); return a; }
    int     prf_gilow1  (void *handle, int a)     { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT,       MPI_MIN, rcfd->Parallel->world); return a; }
    double  prf_grsum1  (void *handle, double a)  { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_DOUBLE,    MPI_SUM, rcfd->Parallel->world); return a; }
    double  prf_grhigh1 (void *handle, double a)  { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_DOUBLE,    MPI_MAX, rcfd->Parallel->world); return a; }
    double  prf_grlow1  (void *handle, double a)  { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_DOUBLE,    MPI_MIN, rcfd->Parallel->world); return a; }
    int64_t prf_gi64sum1(void *handle, int64_t a) { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT64_T,   MPI_SUM, rcfd->Parallel->world); return a; }
    int64_t prf_gi64high1(void *handle,int64_t a) { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT64_T,   MPI_MAX, rcfd->Parallel->world); return a; }
    int64_t prf_gi64low1(void *handle, int64_t a) { rCFD_type *rcfd = (rCFD_type*)handle; MPI_Allreduce(MPI_IN_PLACE, &a, 1, MPI_INT64_T,   MPI_MIN, rcfd->Parallel->world); return a; }
#endif // RCFD_LIB_MPI
#else
    /**
     * @brief Gather varying amounts of integer data from all processes to the root process
     * @param sarray        starting address of send buffer
     * @param scount        number of elements in send buffer
     * @param rarray        address of receive buffer (significant only at root)
     * @param recvcounts    array (of length compute_node_count) containing the number of elements
     *                      that are received from each process (significant only at root)
     * @param displs        array (of length compute_node_count) containing the displacements
     *                      relative to rarray at which to place the incoming data from each process
     *                      (significant only at root)
     * @param root          rank of receiving process
     */
    void prf_gatherv_int(int *sarray, int scount, int *rarray, const int recvcounts[], const int displs[], int root)
    {
#if RP_NODE
        int i_proc;

        if(myid != root){
            if(scount > 0){
                PRF_CSEND_INT(root, sarray, scount, TAG_GATHERV);
            }
            return;
        }

        for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
            if(recvcounts[i_proc] > 0){
                if(i_proc == root){
                    if(scount > 0){
                        memcpy(&rarray[displs[i_proc]], sarray, recvcounts[i_proc]*sizeof(int));
                    }
                }
                else{
                    PRF_CRECV_INT(i_proc, &rarray[displs[i_proc]], recvcounts[i_proc], TAG_GATHERV);
                }
            }
        }
#endif // RP_NODE
    }

    /**
     * @brief Gather varying amounts of floating point data from all processes to the root process
     * @param sarray        starting address of send buffer
     * @param scount        number of elements in send buffer
     * @param rarray        address of receive buffer (significant only at root)
     * @param recvcounts    array (of length compute_node_count) containing the number of elements
     *                      that are received from each process (significant only at root)
     * @param displs        array (of length compute_node_count) containing the displacements
     *                      relative to rarray at which to place the incoming data from each process
     *                      (significant only at root)
     * @param root          rank of receiving process
     */
    void prf_gatherv_real(double *sarray, int scount, double *rarray, const int recvcounts[], const int displs[], int root)
    {
#if RP_NODE
        int i_proc;

        if(myid != root){
            if(scount > 0){
                PRF_CSEND_REAL(root, sarray, scount, TAG_GATHERV);
            }
            return;
        }

        for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
            if(recvcounts[i_proc] > 0){
                if(i_proc == root){
                    if(scount > 0){
                        memcpy(&rarray[displs[i_proc]], sarray, recvcounts[i_proc] * sizeof(double));
                    }
                }
                else{
                    PRF_CRECV_REAL(i_proc, &rarray[displs[i_proc]], recvcounts[i_proc], TAG_GATHERV);
                }
            }
        }
#endif // RP_NODE
    }

#endif // RCFD_STANDALONE

    /**
     * @brief Gather varying amounts of integer data from all processes to process 0,
     *          size of receive buffer is calculated and memory is allocated
     * @param handle    pointer to rCFD instance
     * @param sarray    starting address of send buffer
     * @param scount    number of elements in send buffer
     * @param rarray0   address of receive buffer to be allocated (significant only at root)
     * @return total number of elements received on process 0
     */
    int prf_create_gatherv0_int(void *handle, int *sarray, int scount, int **rarray0)
    {
#if defined(RCFD_STANDALONE) && defined(RCFD_LIB_MPI)
        rCFD_type *rcfd = (rCFD_type*)handle;
#else
        (void)(handle);
#endif
        int count0 = 0;

#if RP_NODE

        int *recvcounts = NULL, *displs = NULL;

        recvcounts = (int*)malloc(sizeof(int) * compute_node_count);

        // gather recvcount for each processor
        PRF_GATHER_INT(recvcounts, &scount, 1, 0);

        if(myid == 0){
            int i_proc;
            // create offsets in recv buffer and compute total recv count
            displs = (int*)malloc(sizeof(int) * compute_node_count);
            displs[0] = 0;
            for(i_proc = 1; i_proc < compute_node_count; ++i_proc){
                count0 += recvcounts[i_proc-1];
                displs[i_proc] = displs[i_proc-1] + recvcounts[i_proc-1];
            }
            count0 += recvcounts[compute_node_count-1];
            // allocate recv buffer
            *rarray0 = (int*)malloc(sizeof(int) * count0);
        }

        PRF_BCAST_INT(&count0, 1, 0);

        PRF_GATHERV_INT(*rarray0, sarray, recvcounts, scount, displs, 0);

        free(recvcounts);
        free(displs);

#endif // RP_NODE

        return count0;
    }

    /**
     * @brief Gather varying amounts of floating point data from all processes to process 0,
     *          size of receive buffer is calculated and memory is allocated
     * @param handle    pointer to rCFD instance
     * @param sarray    starting address of send buffer
     * @param scount    number of elements in send buffer
     * @param rarray0   address of receive buffer to be allocated (significant only at root)
     * @return total number of elements received on process 0
     */
    int prf_create_gatherv0_real(void *handle, double *sarray, int scount, double **rarray0)
    {
#if defined(RCFD_STANDALONE) && defined(RCFD_LIB_MPI)
        rCFD_type *rcfd = (rCFD_type*)handle;
#else
        (void)(handle);
#endif
        int count0 = 0;

#if RP_NODE

        int *recvcounts = NULL, *displs = NULL;

        recvcounts = (int*)malloc(sizeof(int) * compute_node_count);

        // gather recvcount for each processor
        PRF_GATHER_INT(recvcounts, &scount, 1, 0);

        if(myid == 0){
            int i_proc;
            // create offsets in recv buffer and compute total recv count
            displs = (int*)malloc(sizeof(int) * compute_node_count);
            displs[0] = 0;
            for(i_proc = 1; i_proc < compute_node_count; ++i_proc){
                count0 += recvcounts[i_proc-1];
                displs[i_proc] = displs[i_proc-1] + recvcounts[i_proc-1];
            }
            count0 += recvcounts[compute_node_count-1];
            // allocate recv buffer
            *rarray0 = (double*)malloc(sizeof(double) * count0);
        }

        PRF_BCAST_INT(&count0, 1, 0);

        PRF_GATHERV_REAL(*rarray0, sarray, recvcounts, scount, displs, 0);

        free(recvcounts);
        free(displs);

#endif // RP_NODE

        return count0;
    }


    void parallel_init(rCFD_type *rcfd, MPI_Comm communicator)
    {
#if RP_HOST
        rcfd->Parallel->MPI_Topo    = NULL;
        rcfd->Parallel->C2Cs_MPI    = NULL;
#endif
#if RP_NODE
        rcfd->Parallel->MPI_Topo    = (MPI_topo_type*)malloc(sizeof(MPI_topo_type));
        rcfd->Parallel->C2Cs_MPI    = (C2C_MPI_type*)malloc(sizeof(C2C_MPI_type));
        rcfd->Parallel->MPI_Topo->MPI_Cell = NULL;
        rcfd->Parallel->C2Cs_MPI->shifts_in = NULL;
#endif
        rcfd->Parallel->finalize = 0;
        rcfd->Parallel->world = communicator;
        MPI_Comm_rank(communicator, &rcfd->Parallel->rank);
        MPI_Comm_size(communicator, &rcfd->Parallel->size);
    }

#if 1 /* defaults for ANSYS/Fluent */
#if RP_NODE

    void parallel_read_hosts(rCFD_type *rcfd, short i_layer)
    {
        FILE    *f_topo = NULL;

        if(i_layer > 0){

            return;
        }

        {
            char    filename[80];

            snprintf(filename, 80, "%s_%d", rcfd->File_Dict->Topo_filename, myid);

            f_topo = fopen(filename, "r");

            if(f_topo == NULL){

                report_info(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                            "[%d] parallel_read_hosts: failed to open '%s'. "
                            "Expecting topology to be set via user functions",
                            myid, filename);
                return;
            }
        }

        // read cell data
        {
            int     i_cell;
            int     number_of_cells = 0;
            int     number_of_faces = 0;
            int     number_of_nodes = 0;
            double  x[3];
            double  volume;

            fscanf(f_topo, "%d %d %d \n", &number_of_cells, &number_of_faces, &number_of_nodes);

            if(_CELL_DICT(i_layer).number_of_cells > 0 && _CELL_DICT(i_layer).number_of_cells != number_of_cells){

                report_info(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                            "[%d] parallel_read_hosts: mismatch of cell count in Cell_Dict (%d) and Topo file (%d)."
                            "Expecting topology to be set via user functions",
                            myid, _CELL_DICT(i_layer).number_of_cells, number_of_cells);

                fclose(f_topo);

                return;
            }

            //loop_cells
            for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                fscanf(f_topo, "%le %le %le %le %d \n", &x[0], &x[1], &x[2], &volume, &_MPI_Cell.host_of_cell[i_cell]);
            }
        }
    }

    void parallel_set_host_of_cell(rCFD_type *rcfd, short i_layer)
    {
        parallel_read_hosts(rcfd, i_layer);

        rcfd->rCFD_UDF->_user_init_cell_hosts(rcfd, i_layer);
    }

#endif
#endif

#if 1 /* communication tools */
#if RP_NODE

    enum{ /* corona_communication_formats */

        map_host_format,
        sum_cells_format
    };

    void parallel_communicate(rCFD_type *rcfd, double *cell_data, const short i_layer, const short comm_format)
    {
        int     i_cell, i_list, i_node;

        int     size = 0;

        double  *list_of_double;

        /* MPI.1 Node-0 fills _MPI_Cell.data with corona cell values */
        {
            if(myid > 0){

                list_of_double = (double*)malloc(_CELL_DICT(i_layer).number_of_ext_cells * sizeof(double));

                i_list = 0;

                //loop_ext_cells
                for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    list_of_double[i_list] = cell_data[i_cell];

                    i_list++;
                }

                PRF_CSEND_INT(node_zero, &_CELL_DICT(i_layer).number_of_ext_cells, 1, myid);

                PRF_CSEND_REAL(node_zero, list_of_double, _CELL_DICT(i_layer).number_of_ext_cells, myid);

                free(list_of_double);
            }

            if(myid == 0){

                i_list = 0;

                //loop_ext_cells
                for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    _MPI_Cell.data[i_list] = cell_data[i_cell];

                    i_list++;
                }

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CRECV_INT(i_node, &size, 1, i_node);

                    list_of_double = (double*)malloc(size * sizeof(double));

                    PRF_CRECV_REAL(i_node, list_of_double, size, i_node);

                    for(i_cell = 0; i_cell < size; i_cell++){

                        _MPI_Cell.data[i_list] = list_of_double[i_cell];

                        i_list++;
                    }

                    free(list_of_double);
                }
            }
        }

        /* MPI.2 Node-0 sends _MPI_Cell.data to host-Nodes */
        {
            if(myid == 0){

                i_list = _MPI_Cell.number_of_host_cells_per_node[0];

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CSEND_INT(i_node, &_MPI_Cell.number_of_host_cells_per_node[i_node], 1, node_zero);

                    list_of_double = (double*)malloc(_MPI_Cell.number_of_host_cells_per_node[i_node] * sizeof(double));

                    for(i_cell = 0; i_cell < _MPI_Cell.number_of_host_cells_per_node[i_node]; i_cell++){

                        list_of_double[i_cell] = _MPI_Cell.data[_MPI_Cell.host2ext_index[i_list]];

                        i_list++;
                    }

                    PRF_CSEND_REAL(i_node, list_of_double, _MPI_Cell.number_of_host_cells_per_node[i_node], node_zero);

                    free(list_of_double);
                }
            }

        }

        /* MPI.3 all nodes update data based on hosting inner cells; Node-0 collects updated data into _MPI_Cell.data */
        {

            /* other nodes receive corona cells, eventually update their host cells and send them back to node-0 */
            if(myid > 0){

                PRF_CRECV_INT(node_zero, &size, 1, node_zero);

                list_of_double = (double*)malloc(size * sizeof(double));

                PRF_CRECV_REAL(node_zero, list_of_double, size, node_zero);

                // eventually update host_cell data by corona cell values
                if(comm_format == sum_cells_format){

                    for(i_list = 0; i_list < size; i_list++){

                        cell_data[_MPI_Cell.hosting_cell_index[i_list]] += list_of_double[i_list];
                    }
                }

                // send host cell values back to node-0
                for(i_list = 0; i_list < size; i_list++){

                    list_of_double[i_list] = cell_data[_MPI_Cell.hosting_cell_index[i_list]];
                }

                PRF_CSEND_INT(node_zero, &size, 1, myid);

                PRF_CSEND_REAL(node_zero, list_of_double, size, myid);

                free(list_of_double);
            }


            if(myid == 0){

                /* Node-0 fills his own _MPI_Cell.data */

                size = _MPI_Cell.number_of_host_cells_per_node[0];

                list_of_double = (double*)malloc(size * sizeof(double));

                // fill list of host cells of node-0
                for(i_cell = 0; i_cell < size; i_cell++){

                    list_of_double[i_cell] = _MPI_Cell.data[_MPI_Cell.host2ext_index[i_cell]];
                }

                //
                if(comm_format == sum_cells_format){

                    for(i_list = 0; i_list < size; i_list++){

                        cell_data[_MPI_Cell.hosting_cell_index[i_list]] += list_of_double[i_list];
                    }
                }

                for(i_list = 0; i_list < size; i_list++){

                    list_of_double[i_list] = cell_data[_MPI_Cell.hosting_cell_index[i_list]];
                }

                for(i_cell = 0; i_cell < size; i_cell++){

                    _MPI_Cell.data[_MPI_Cell.host2ext_index[i_cell]] = list_of_double[i_cell];
                }

                free(list_of_double);

                /* Node-0 collects _MPI_Cell.data from other nodes */

                i_list = _MPI_Cell.number_of_host_cells_per_node[0];

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CRECV_INT(i_node, &size, 1, i_node);

                    list_of_double = (double*)malloc(size * sizeof(double));

                    PRF_CRECV_REAL(i_node, list_of_double, size, i_node);

                    for(i_cell = 0; i_cell < size; i_cell++){

                        _MPI_Cell.data[_MPI_Cell.host2ext_index[i_list]] = list_of_double[i_cell];

                        i_list++;
                    }

                    free(list_of_double);
                }
            }

        }

        /* MPI.4 Node-0 sends updated data to nodes */
        {

            if(myid == 0){

                i_list = 0;

                //loop_ext_cells
                for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    cell_data[i_cell] = _MPI_Cell.data[i_list];

                    i_list++;
                }

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CSEND_INT(i_node, &_MPI_Cell.number_of_ext_cells_per_node[i_node], 1, node_zero);

                    list_of_double = (double*)malloc(_MPI_Cell.number_of_ext_cells_per_node[i_node] * sizeof(double));

                    for(i_cell = 0; i_cell < _MPI_Cell.number_of_ext_cells_per_node[i_node]; i_cell++){

                        list_of_double[i_cell] = _MPI_Cell.data[i_list];

                        i_list++;
                    }

                    PRF_CSEND_REAL(i_node, list_of_double, _MPI_Cell.number_of_ext_cells_per_node[i_node], node_zero);

                    free(list_of_double);
                }

            }

            if(myid > 0){

                PRF_CRECV_INT(node_zero, &size, 1, node_zero);

                list_of_double = (double*)malloc(size * sizeof(double));

                PRF_CRECV_REAL(node_zero, list_of_double, size, node_zero);

                i_list = 0;

                //loop_ext_cells
                for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                    cell_data[i_cell] = list_of_double[i_list];

                    i_list++;
                }

                free(list_of_double);
            }
        }

    }

    void parallel_communicate_data(rCFD_type *rcfd, const short i_layer, const short i_phase, const short i_data, const short comm_format)
    {
        int     i_cell;

        // fill _.tmp field with host and corona values
        //loop_cells
        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

            _C(i_layer).tmp[i_cell] = _C(i_layer).data[i_phase][i_cell][i_data];
        }

        parallel_communicate(rcfd, _C(i_layer).tmp, i_layer, comm_format);

        // re-fill ext cells with updated corona values
        //loop_ext_cells
        for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

            _C(i_layer).data[i_phase][i_cell][i_data] = _C(i_layer).tmp[i_cell];
        }
    }

#endif
#endif

#if 1 /* init */
#if RP_NODE

    void parallel_init_topology(rCFD_type *rcfd, short i_layer)
    {

#if 1   // local vars
        int         i_node, i_cell, i_list, i_host;

        int         size;

        short       local_debug_on = 0;

        int         *list_of_int = NULL;

        double      *list_of_double = NULL;
#endif

        // only implemented for base layer !
        if(i_layer > 0){

            _MPI_Cell.total_number_of_ext_cells = 0;

            _MPI_Cell.number_of_ext_cells_per_node = NULL;

            _MPI_Cell.number_of_host_cells_per_node = NULL;

            _MPI_Cell.host_of_cell = NULL;

            _MPI_Cell.host_of_ext_cells = NULL;

            _MPI_Cell.hosting_cell_index = NULL;

            _MPI_Cell.host2ext_index = NULL;

            _MPI_Cell.data = NULL;

            return;
        }

        /* P.1  allocate MPI_topo and _MPI_Cell  */
        {

            rcfd->Parallel->MPI_Topo->MPI_Cell = (MPI_cells_type*)malloc(rcfd->Solver_Dict->number_of_layers * sizeof(MPI_cells_type));

            _MPI_Cell.host_of_cell = (int*)malloc(_CELL_DICT(i_layer).number_of_cells * sizeof(int));

            _MPI_Cell.hosting_cell_index = NULL;

            if(myid > 0){

                PRF_CSEND_INT(node_zero, &_CELL_DICT(i_layer).number_of_ext_cells, 1, myid);

                _MPI_Cell.total_number_of_ext_cells = -1;

                _MPI_Cell.number_of_ext_cells_per_node = NULL;

                _MPI_Cell.number_of_host_cells_per_node = NULL;

                _MPI_Cell.host_of_ext_cells = NULL;

                _MPI_Cell.host2ext_index = NULL;

                _MPI_Cell.data = NULL;
            }

            if(myid == 0){

                // per node and total number corona cells

                _MPI_Cell.number_of_ext_cells_per_node = (int*)malloc(compute_node_count * sizeof(int));

                _MPI_Cell.number_of_ext_cells_per_node[0] = _CELL_DICT(i_layer).number_of_ext_cells;

                _MPI_Cell.total_number_of_ext_cells = _CELL_DICT(i_layer).number_of_ext_cells;

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CRECV_INT(i_node, &_MPI_Cell.number_of_ext_cells_per_node[i_node], 1, i_node);

                    _MPI_Cell.total_number_of_ext_cells += _MPI_Cell.number_of_ext_cells_per_node[i_node];
                }

                // other node-0 MPI vars allocated, but not set yet

                _MPI_Cell.number_of_host_cells_per_node = (int*)malloc(compute_node_count * sizeof(int));

                _MPI_Cell.host_of_ext_cells = (int*)malloc(_MPI_Cell.total_number_of_ext_cells * sizeof(int));

                _MPI_Cell.host2ext_index = (int*)malloc(_MPI_Cell.total_number_of_ext_cells * sizeof(int));

                _MPI_Cell.data = (double*)malloc(_MPI_Cell.total_number_of_ext_cells * sizeof(double));
            }
        }

        /* P.2  all nodes set _MPI_Cell.host_of_cell */
        {
            parallel_set_host_of_cell(rcfd, i_layer);
        }

        /* P.3  set up MPI_cell communication */
        {
            double *tmp_x;

            /* P3.1 node-0 fills _MPI_Cell.host_of_ext_cells, tmp_x */
            {
                if(myid > 0){

                    list_of_int = (int*)malloc( _CELL_DICT(i_layer).number_of_ext_cells * sizeof(int));

                    list_of_double = (double*)malloc( 3 * _CELL_DICT(i_layer).number_of_ext_cells * sizeof(double));

                    i_list = 0;

                    //loop_ext_cells
                    for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                        list_of_int[i_list] = _MPI_Cell.host_of_cell[i_cell];

                        list_of_double[(3 * i_list)] =      _C(i_layer).x[i_cell][0];
                        list_of_double[(3 * i_list + 1)] =  _C(i_layer).x[i_cell][1];
                        list_of_double[(3 * i_list + 2)] =  _C(i_layer).x[i_cell][2];

                        i_list++;
                    }

                    PRF_CSEND_INT(node_zero, list_of_int, _CELL_DICT(i_layer).number_of_ext_cells, myid);

                    PRF_CSEND_REAL(node_zero, list_of_double, (3 * _CELL_DICT(i_layer).number_of_ext_cells), myid);

                    free(list_of_int);

                    free(list_of_double);
                }

                if(myid == 0){

                    tmp_x = (double*)malloc(3 * _MPI_Cell.total_number_of_ext_cells * sizeof(double));

                    i_list = 0;

                    //loop_ext_cells
                    for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                        _MPI_Cell.host_of_ext_cells[i_list] = _MPI_Cell.host_of_cell[i_cell];

                        tmp_x[(3 * i_list)] =       _C(i_layer).x[i_cell][0];
                        tmp_x[(3 * i_list + 1)] =   _C(i_layer).x[i_cell][1];
                        tmp_x[(3 * i_list + 2)] =   _C(i_layer).x[i_cell][2];

                        i_list++;
                    }

                    /* get ext cell info from Nodes */
                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        list_of_int = (int*)malloc(_MPI_Cell.number_of_ext_cells_per_node[i_node] * sizeof(int));

                        list_of_double = (double*)malloc( 3 * _MPI_Cell.number_of_ext_cells_per_node[i_node] * sizeof(double));

                        PRF_CRECV_INT(i_node, list_of_int, _MPI_Cell.number_of_ext_cells_per_node[i_node], i_node);

                        PRF_CRECV_REAL(i_node, list_of_double, (3 * _MPI_Cell.number_of_ext_cells_per_node[i_node]), i_node);

                        for(i_cell = 0; i_cell < _MPI_Cell.number_of_ext_cells_per_node[i_node]; i_cell++){

                            _MPI_Cell.host_of_ext_cells[i_list] = list_of_int[i_cell];

                            tmp_x[(3 * i_list)] =       list_of_double[(3 * i_cell)];
                            tmp_x[(3 * i_list + 1)] =   list_of_double[(3 * i_cell + 1)];
                            tmp_x[(3 * i_list + 2)] =   list_of_double[(3 * i_cell + 2)];

                            i_list++;
                        }

                        free(list_of_int);

                        free(list_of_double);
                    }
                }
            }

            /* P3.2 analyze _MPI_Cell.host_of_ext_cells, set _MPI_Cell.host2ext_index */
            {
                if(myid == 0){

                    for(i_host = 0; i_host < compute_node_count; i_host++){

                        _MPI_Cell.number_of_host_cells_per_node[i_host] = 0;
                    }

                    // set number of host cells per node
                    for(i_cell = 0; i_cell < _MPI_Cell.total_number_of_ext_cells; i_cell++){

                        i_host = _MPI_Cell.host_of_ext_cells[i_cell];

                        _MPI_Cell.number_of_host_cells_per_node[i_host]++;
                    }

                    int host_start_index[compute_node_count], host_index[compute_node_count];

                    host_start_index[0] = 0; host_index[0] = 0;

                    for(i_host = 1; i_host < compute_node_count; i_host++){

                        host_start_index[i_host] = host_start_index[(i_host-1)] + _MPI_Cell.number_of_host_cells_per_node[(i_host-1)];

                        host_index[i_host] = 0;
                    }

                    for(i_cell = 0; i_cell < _MPI_Cell.total_number_of_ext_cells; i_cell++){

                        i_host = _MPI_Cell.host_of_ext_cells[i_cell];

                        _MPI_Cell.host2ext_index[(host_start_index[i_host] + host_index[i_host])] = i_cell;

                        host_index[i_host]++;
                    }
                }
            }

            /* P3.3 set _MPI_Cell.hosting_cell_index */
            {
                if(myid == 0){

                    _MPI_Cell.number_of_host_cells = _MPI_Cell.number_of_host_cells_per_node[0];

                    list_of_double = (double*)malloc( 3 * _MPI_Cell.number_of_host_cells_per_node[0] * sizeof(double));

                    i_list = 0;

                    for(i_cell = 0; i_cell < _MPI_Cell.total_number_of_ext_cells; i_cell++){

                        i_host = _MPI_Cell.host_of_ext_cells[i_cell];

                        if(i_host == 0){

                            list_of_double[(3 * i_list)] =      tmp_x[(3 * i_cell)];
                            list_of_double[(3 * i_list + 1)] =  tmp_x[(3 * i_cell + 1)];
                            list_of_double[(3 * i_list + 2)] =  tmp_x[(3 * i_cell + 2)];

                            i_list++;
                        }
                    }

                    _MPI_Cell.hosting_cell_index = (int*)malloc(_MPI_Cell.number_of_host_cells * sizeof(int));

                    for(i_host = 0; i_host < _MPI_Cell.number_of_host_cells; i_host++){

                        _MPI_Cell.hosting_cell_index[i_host] = -1;
                    }

                    for(i_list = 0; i_list < _MPI_Cell.number_of_host_cells_per_node[0]; i_list++){

                        //loop_int_cells
                        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                            //i_cell = _F(i_layer).c0[i_face];     // potential host cell

                            if( (_C(i_layer).x[i_cell][0] == list_of_double[(3 * i_list + 0)]) &&
                                (_C(i_layer).x[i_cell][1] == list_of_double[(3 * i_list + 1)]) &&
                                (_C(i_layer).x[i_cell][2] == list_of_double[(3 * i_list + 2)])){

                                _MPI_Cell.hosting_cell_index[i_list] = i_cell;
                            }
                        }
                    }

                    free(list_of_double);

                    /* send tmp_x to Nodes */
                    for(i_node = 1; i_node < compute_node_count; i_node++){

                        PRF_CSEND_INT(i_node, &_MPI_Cell.number_of_host_cells_per_node[i_node], 1, node_zero);

                        list_of_double = (double*)malloc( 3 * _MPI_Cell.number_of_host_cells_per_node[i_node] * sizeof(double));

                        i_list = 0;

                        for(i_cell = 0; i_cell < _MPI_Cell.total_number_of_ext_cells; i_cell++){

                            i_host = _MPI_Cell.host_of_ext_cells[i_cell];

                            if(i_host == i_node){

                                list_of_double[(3 * i_list)] =      tmp_x[(3 * i_cell)];
                                list_of_double[(3 * i_list + 1)] =  tmp_x[(3 * i_cell + 1)];
                                list_of_double[(3 * i_list + 2)] =  tmp_x[(3 * i_cell + 2)];

                                i_list++;
                            }
                        }

                        PRF_CSEND_REAL(i_node, list_of_double, (3 * _MPI_Cell.number_of_host_cells_per_node[i_node]), node_zero);

                        free(list_of_double);
                    }

                    free(tmp_x);
                }

                if(myid > 0){

                    PRF_CRECV_INT(node_zero, &size, 1, node_zero);  /* size = number of host cells for this node */

                    _MPI_Cell.number_of_host_cells = size;

                    _MPI_Cell.hosting_cell_index = (int*)malloc(size * sizeof(int));

                    for(i_host = 0; i_host < _MPI_Cell.number_of_host_cells; i_host++){

                        _MPI_Cell.hosting_cell_index[i_host] = -1;
                    }

                    list_of_double = (double*)malloc(3 * size * sizeof(double));

                    PRF_CRECV_REAL(node_zero, list_of_double, 3 * size, node_zero);

                    for(i_list = 0; i_list < size; i_list++){

                        //loop_int_cells
                        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_int_cells; i_cell++){

                            //i_cell = _F(i_layer).c0[i_face];     // potential host cell

                            if( (_C(i_layer).x[i_cell][0] == list_of_double[(3 * i_list + 0)]) &&
                                (_C(i_layer).x[i_cell][1] == list_of_double[(3 * i_list + 1)]) &&
                                (_C(i_layer).x[i_cell][2] == list_of_double[(3 * i_list + 2)])){

                                _MPI_Cell.hosting_cell_index[i_list] = i_cell;
                            }
                        }
                    }

                    free(list_of_double);
                }

                if(local_debug_on){

                    for(i_host = 0; i_host < _MPI_Cell.number_of_host_cells; i_host++){

                        if((_MPI_Cell.hosting_cell_index[i_host] < 0)||(_MPI_Cell.hosting_cell_index[i_host] > _CELL_DICT(i_layer).number_of_int_cells)){

                            report_debug(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                                         "[%d] _MPI_Cell.hosting_cell_index[%d] = %d",
                                         myid, i_host, _MPI_Cell.hosting_cell_index[i_host]);
                        }
                    }
                }
            }

        }

        if(local_debug_on){

            if(myid == 0){

                report_debug0(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                              "total_number_of_ext_cells = %d",
                              _MPI_Cell.total_number_of_ext_cells);

                for(i_node = 0; i_node < compute_node_count; i_node++){

                    report_debug0(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                                  "number_of_ext/host_cells_per_node[%d] = %d / %d",
                                  i_node,
                                  _MPI_Cell.number_of_ext_cells_per_node[i_node],
                                  _MPI_Cell.number_of_host_cells_per_node[i_node]);
                }
            }
        }

    }

    void parallel_expand_initial_fields(rCFD_type *rcfd, short i_layer)
    {

        int     i_rec_user, i_frame, i_cell;

        double  *tmp = NULL;

        if(_C(i_layer).rec_user != NULL){

            // expand rec_user values to corona cells
            {
                tmp = (double*)malloc(_CELL_DICT(i_layer).number_of_cells * sizeof(double));

                //loop_frames
                for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                    for(i_rec_user = 0; i_rec_user < _CELL_DICT(i_layer).number_of_rec_user_vars; i_rec_user++){

                        //loop_cells
                        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                            tmp[i_cell] = _C(i_layer).rec_user[i_frame][i_cell][i_rec_user];
                        }

                        parallel_communicate(rcfd, tmp, i_layer, map_host_format);

                        //loop_ext_cells
                        for(i_cell = _CELL_DICT(i_layer).number_of_int_cells; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                            _C(i_layer).rec_user[i_frame][i_cell][i_rec_user] = tmp[i_cell];
                        }
                    }
                }

                free(tmp);
            }
        }
    }

    void init_parallel_C2Cs(rCFD_type *rcfd)
    {

        int                 i_state, i_phase, i_frame, i_layer, i_out, i_node, i_shift;
        int                 number_of_shifts, max_number_of_shifts, index_in_total_list;
        int                 *list_of_int = NULL;
        double              *list_of_double = NULL;

        tmp_C2C_type        tmp_C2Cs;

        /* A. allocate C2C_MPI.shifts (shifts.data is allocated later) */
        {

            max_number_of_shifts = 0;

            //loop_states
            for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                //loop_phases
                for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                    //loop_frames
                    for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                        //loop_layers
                        for(i_layer = 0; i_layer < rcfd->Solver_Dict->number_of_layers; i_layer++){

                            number_of_shifts = PRF_GISUM1(rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer]);

                            if(number_of_shifts > max_number_of_shifts){

                                max_number_of_shifts = number_of_shifts;
                            }
                        }
                    }
                }
            }

            if(myid == 0){

                rcfd->Parallel->C2Cs_MPI->shifts_in = (C2C_MPI_shift_type*)malloc(max_number_of_shifts * sizeof(C2C_MPI_shift_type));
            }
            else{
                rcfd->Parallel->C2Cs_MPI->shifts_in = NULL;
            }

            rcfd->Parallel->C2Cs_MPI->max_number_of_MPI_shifts = max_number_of_shifts;

            rcfd->Parallel->C2Cs_MPI->number_of_shifts_in = 0;
        }

        /* B. determine C2Cs.shift_out, C2Cs.in2out, C2Cs.number_of_shifts_to_node_zero */
        {

            //loop_states
            for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                //loop_phases
                for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                    //loop_frames
                    for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                        rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_to_node_zero =    (int**)malloc(rcfd->Solver_Dict->number_of_layers * sizeof(int*));
                        rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_from_node_zero =  (int**)malloc(rcfd->Solver_Dict->number_of_layers * sizeof(int*));
                        rcfd->C2Cs[i_state][i_phase][i_frame].in2out =                           (int***)malloc(rcfd->Solver_Dict->number_of_layers * sizeof(int**));

                        //loop_layers
                        for(i_layer = 0; i_layer < rcfd->Solver_Dict->number_of_layers; i_layer++){

                            /* B.1: send local rCFD_C2C_incoming entries to Node-0,      */
                            /*      which collects them into rCFD_C2C_MPI_buffer         */
                            {

                                if (myid > 0){

                                    number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];

                                    PRF_CSEND_INT(node_zero, &number_of_shifts, 1, myid);

                                    if(number_of_shifts > 0){

                                        list_of_int=(int*)malloc(number_of_shifts * sizeof(int));
                                        list_of_double=(double*)malloc(number_of_shifts * sizeof(double));

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            list_of_int[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;
                                        }

                                        PRF_CSEND_INT(node_zero, list_of_int, number_of_shifts, myid);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            list_of_int[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;
                                        }

                                        PRF_CSEND_INT(node_zero, list_of_int, number_of_shifts, myid);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            list_of_int[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;
                                        }

                                        PRF_CSEND_INT(node_zero, list_of_int, number_of_shifts, myid);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            list_of_int[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node1;
                                        }

                                        PRF_CSEND_INT(node_zero, list_of_int, number_of_shifts, myid);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            list_of_double[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;
                                        }

                                        PRF_CSEND_REAL(node_zero, list_of_double, number_of_shifts, myid);

                                        free(list_of_int);
                                        free(list_of_double);
                                    }
                                }

                                if (myid == 0){

                                    /* initiate for each frame */
                                    for(i_shift = 0; i_shift < rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer]; i_shift++){

                                        rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c0     = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;
                                        rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0  = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;
                                        rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c1     = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;
                                        rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node1  = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node1;
                                        rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].w0     = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;
                                    }

                                    index_in_total_list = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];    /* first index from Node-0 */

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_to_node_zero[i_layer] = (int*)malloc(compute_node_count * sizeof(int));

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_to_node_zero[i_layer][0] = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];

                                    for(i_node = 1; i_node < compute_node_count; i_node++){

                                        PRF_CRECV_INT(i_node, &number_of_shifts, 1, i_node);

                                        if(number_of_shifts > 0){

                                            list_of_int=(int*)malloc(number_of_shifts * sizeof(int));
                                            list_of_double=(double*)malloc(number_of_shifts * sizeof(double));

                                            PRF_CRECV_INT(i_node, list_of_int, number_of_shifts, i_node);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                rcfd->Parallel->C2Cs_MPI->shifts_in[index_in_total_list + i_shift].c0 = list_of_int[i_shift];
                                            }

                                            PRF_CRECV_INT(i_node, list_of_int, number_of_shifts, i_node);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                rcfd->Parallel->C2Cs_MPI->shifts_in[index_in_total_list + i_shift].node0 = list_of_int[i_shift];
                                            }

                                            PRF_CRECV_INT(i_node, list_of_int, number_of_shifts, i_node);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                rcfd->Parallel->C2Cs_MPI->shifts_in[index_in_total_list + i_shift].c1 = list_of_int[i_shift];
                                            }

                                            PRF_CRECV_INT(i_node, list_of_int, number_of_shifts, i_node);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                rcfd->Parallel->C2Cs_MPI->shifts_in[index_in_total_list + i_shift].node1 = list_of_int[i_shift];
                                            }

                                            PRF_CRECV_REAL(i_node, list_of_double, number_of_shifts, i_node);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                rcfd->Parallel->C2Cs_MPI->shifts_in[index_in_total_list + i_shift].w0 = list_of_double[i_shift];
                                            }

                                            free(list_of_int);
                                            free(list_of_double);
                                        }

                                        rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_to_node_zero[i_layer][i_node] = number_of_shifts;

                                        index_in_total_list += number_of_shifts;
                                    }

                                    rcfd->Parallel->C2Cs_MPI->number_of_shifts_in = index_in_total_list;
                                }
                            }

                            /* B.2: Node-0 defines and sends C2Cs.shifts_out for each Node */
                            {

                                if (myid == 0){

                                    i_out = 0;

                                    for(i_shift = 0; i_shift < rcfd->Parallel->C2Cs_MPI->number_of_shifts_in; i_shift++){

                                        if(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0 == node_zero){

                                            i_out++;
                                        }
                                    }

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_out[i_layer] = i_out;

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer] = (C2C_shift_type*)malloc(i_out * sizeof(C2C_shift_type));

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_from_node_zero[i_layer] = (int*)malloc(compute_node_count * sizeof(int));

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_from_node_zero[i_layer][0] = i_out;


                                    i_out = 0;

                                    /* fill C2C.shifts_out */
                                    for(i_shift = 0; i_shift < rcfd->Parallel->C2Cs_MPI->number_of_shifts_in; i_shift++){

                                        if(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0 == node_zero){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_out].c0    = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_out].node0 = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_out].c1    = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_out].node1 = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node1;
                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_out].w0    = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].w0;

                                            i_out++;
                                        }
                                    }

                                    /* loop other nodes */
                                    for(i_node = 1; i_node < compute_node_count; i_node++){

                                        i_out = 0;

                                        /* get number of outgoing tracers */
                                        for(i_shift = 0; i_shift < rcfd->Parallel->C2Cs_MPI->number_of_shifts_in; i_shift++){

                                            if(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0 == i_node){

                                                i_out++;
                                            }
                                        }

                                        rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_from_node_zero[i_layer][i_node] = i_out;

                                        /* wrong !! C2Cs[i_state][i_phase][i_frame].number_of_shifts_out = i_out;*/

                                        number_of_shifts = i_out;

                                        /* send outgoing tracers to individual nodes */
                                        PRF_CSEND_INT(i_node, &number_of_shifts, 1, node_zero);

                                        if (number_of_shifts > 0){

                                            /* allocate tmp storage for outgoing tracers */
                                            tmp_C2Cs.shifts_out = (C2C_shift_type*)malloc(number_of_shifts * sizeof(C2C_shift_type));

                                            /* fill tmp storage for outgoing tracers */
                                            i_out = 0;

                                            for(i_shift = 0; i_shift < rcfd->Parallel->C2Cs_MPI->number_of_shifts_in; i_shift++){

                                                if(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0 == i_node){

                                                    tmp_C2Cs.shifts_out[i_out].c0       = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c0;
                                                    tmp_C2Cs.shifts_out[i_out].node0    = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0;
                                                    tmp_C2Cs.shifts_out[i_out].c1       = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].c1;
                                                    tmp_C2Cs.shifts_out[i_out].node1    = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node1;
                                                    tmp_C2Cs.shifts_out[i_out].w0       = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].w0;

                                                    i_out++;
                                                }
                                            }

                                            list_of_int=(int*)malloc(number_of_shifts * sizeof(int));
                                            list_of_double=(double*)malloc(number_of_shifts * sizeof(double));

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                list_of_int[i_shift] = tmp_C2Cs.shifts_out[i_shift].c0;
                                            }

                                            PRF_CSEND_INT(i_node, list_of_int, number_of_shifts, node_zero);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                list_of_int[i_shift] = tmp_C2Cs.shifts_out[i_shift].node0;
                                            }

                                            PRF_CSEND_INT(i_node, list_of_int, number_of_shifts, node_zero);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                list_of_int[i_shift] = tmp_C2Cs.shifts_out[i_shift].c1;
                                            }

                                            PRF_CSEND_INT(i_node, list_of_int, number_of_shifts, node_zero);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                list_of_int[i_shift] = tmp_C2Cs.shifts_out[i_shift].node1;
                                            }

                                            PRF_CSEND_INT(i_node, list_of_int, number_of_shifts, node_zero);

                                            //loop_shifts
                                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                                list_of_double[i_shift] = tmp_C2Cs.shifts_out[i_shift].w0;
                                            }

                                            PRF_CSEND_REAL(i_node, list_of_double, number_of_shifts, node_zero);

                                            free(list_of_int);
                                            free(list_of_double);

                                            free(tmp_C2Cs.shifts_out);
                                        }

                                    } /* loop other nodes */

                                } /* node == 0 */

                                if (myid > 0){

                                    /* receive rCFD_Tracer_outgoing from Node-0 */

                                    PRF_CRECV_INT(node_zero, &number_of_shifts, 1, node_zero);

                                    rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_out[i_layer] = number_of_shifts;

                                    if(number_of_shifts > 0){

                                        rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer] = (C2C_shift_type*)malloc(number_of_shifts * sizeof(C2C_shift_type));

                                        list_of_int=(int*)malloc(number_of_shifts * sizeof(int));
                                        list_of_double=(double*)malloc(number_of_shifts * sizeof(double));

                                        PRF_CRECV_INT(node_zero, list_of_int, number_of_shifts, node_zero);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].c0 = list_of_int[i_shift];
                                        }

                                        PRF_CRECV_INT(node_zero, list_of_int, number_of_shifts, node_zero);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].node0 = list_of_int[i_shift];
                                        }

                                        PRF_CRECV_INT(node_zero, list_of_int, number_of_shifts, node_zero);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].c1 = list_of_int[i_shift];
                                        }

                                        PRF_CRECV_INT(node_zero, list_of_int, number_of_shifts, node_zero);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].node1 = list_of_int[i_shift];
                                        }

                                        PRF_CRECV_REAL(node_zero, list_of_double, number_of_shifts, node_zero);

                                        //loop_shifts
                                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].w0 = list_of_double[i_shift];
                                        }

                                        free(list_of_int);
                                        free(list_of_double);

                                    }

                                } /* node > 0 */

                            }

                            /* B.3: Node-0 link outgoing data to list of incoming data */
                            {

                                if(myid == 0){

                                    rcfd->C2Cs[i_state][i_phase][i_frame].in2out[i_layer] = (int**)malloc(compute_node_count * sizeof(int*));

                                    for(i_node = 0; i_node < compute_node_count; i_node++){

                                        rcfd->C2Cs[i_state][i_phase][i_frame].in2out[i_layer][i_node] = (int*)malloc(rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_from_node_zero[i_layer][i_node] * sizeof(int));
                                    }


                                    for(i_node = 0; i_node < compute_node_count; i_node++){

                                        i_out = 0;

                                        for(i_shift = 0; i_shift < rcfd->Parallel->C2Cs_MPI->number_of_shifts_in; i_shift++){

                                            if(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].node0 == i_node){

                                                rcfd->C2Cs[i_state][i_phase][i_frame].in2out[i_layer][i_node][i_out] = i_shift;

                                                i_out++;
                                            }
                                        }
                                    }
                                }
                            }

                        } /* loop_layers */

                    } /* loop frames */

                } /* loop phases */

            } /* loop states */

        }

        /* C. Allocate C2C_MPI_shifts.data */
        {
            int     max_number_of_data = 0;

            //loop_phases
            for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                if(max_number_of_data < rcfd->Phase_Dict[i_phase].number_of_data){

                    max_number_of_data = rcfd->Phase_Dict[i_phase].number_of_data;
                }
            }

            if(myid == 0){

                number_of_shifts = rcfd->Parallel->C2Cs_MPI->max_number_of_MPI_shifts;

                //loop_shifts
                for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                    rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].data = (double*)malloc(max_number_of_data * sizeof(double));
                }
            }
        }

    }

    void init_upper_layer_parallel_C2Cs(rCFD_type* rcfd)
    {
        short   debug_this_code = 0;

        if(rcfd->Solver_Dict->number_of_layers > 1){

#if 1       /* local vars */

            int     i_state, i_phase, i_frame, i_layer, i_node, i_cell, i_shift, i_tmp;

            int     number_of_node_cells = 0, number_of_shifts, number_of_shifts_in_upper_layer, number_of_shifts_from_node = 0;

            int     node0, c0, c1, upper_c0, upper_c1;

            int     **tmp_parent_ids = NULL, *list_of_c0, *list_of_node0, *list_of_upper_c0;

            double  w0;
#endif

            i_layer = 0;

            i_tmp = 0;

            while((i_layer + 1) < rcfd->Solver_Dict->number_of_layers){

                /* E.1 Node-0 collects parent_id's */
                {
                    if(myid == 0){

                        tmp_parent_ids = (int**)malloc(compute_node_count * sizeof(int*));

                        tmp_parent_ids[0] = (int*)malloc(_CELL_DICT(i_layer).number_of_cells * sizeof(int));

                        //loop_cells
                        for(i_cell = 0; i_cell < _CELL_DICT(i_layer).number_of_cells; i_cell++){

                            tmp_parent_ids[0][i_cell] = _C(i_layer).parent_cell[i_cell];
                        }

                        for(i_node = 1; i_node < compute_node_count; i_node++){

                            PRF_CRECV_INT(i_node, &number_of_node_cells, 1, i_node);

                            tmp_parent_ids[i_node] = (int*)malloc(number_of_node_cells * sizeof(int));

                            list_of_upper_c0 = (int*)malloc(number_of_node_cells * sizeof(int));

                            PRF_CRECV_INT(i_node, list_of_upper_c0, number_of_node_cells, i_node);

                            for(i_cell = 0; i_cell < number_of_node_cells; i_cell++){

                                tmp_parent_ids[i_node][i_cell] = list_of_upper_c0[i_cell];
                            }

                            free(list_of_upper_c0);
                        }
                    }

                    if(myid > 0){

                        PRF_CSEND_INT(node_zero, &_CELL_DICT(i_layer).number_of_cells, 1, myid);

                        PRF_CSEND_INT(node_zero, _C(i_layer).parent_cell, _CELL_DICT(i_layer).number_of_cells, myid);
                    }

                    if(debug_this_code){

                        report_debug(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE, "[%d] E.1", myid);
                    }
                }

                /* E.2 set (i_layer + 1) C2Cs */
                /* TODO remove doublets */
                //loop_states
                for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                    //loop_phases
                    for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                        //loop_frames
                        for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                            number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];

                            number_of_shifts_in_upper_layer = number_of_shifts;

                            rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[(i_layer + 1)] = number_of_shifts;

                            /* TODO: account for multiple islands */
                            rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[(i_layer + 1)][0] = 0;
                            rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[(i_layer + 1)][1] = number_of_shifts;

                            rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)] = (C2C_shift_type*)malloc(number_of_shifts_in_upper_layer * sizeof(C2C_shift_type));

                            if(myid == 0){

                                /* set (i_layer + 1) shifts_in for node-0 */
                                //loop_shifts
                                for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                    node0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;

                                    c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;

                                    upper_c0 = tmp_parent_ids[node0][c0];

                                    c1 =  rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;

                                    upper_c1 = _C(i_layer).parent_cell[c1];

                                    w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;

                                    if(w0 == 0.0){

                                        i_tmp++;
                                    }

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].node0 = node0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].c0 =    upper_c0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].node1 = myid;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].c1 =    upper_c1;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].w0 =    w0;
                                }

                                /* send list_of_upper_c0 to nodes */
                                for(i_node = 1; i_node < compute_node_count; i_node++){

                                    PRF_CRECV_INT(i_node, &number_of_shifts_from_node, 1, i_node);

                                    list_of_node0 = (int*)malloc(number_of_shifts_from_node * sizeof(int));

                                    list_of_c0 = (int*)malloc(number_of_shifts_from_node * sizeof(int));

                                    PRF_CRECV_INT(i_node, list_of_node0, number_of_shifts_from_node, i_node);

                                    PRF_CRECV_INT(i_node, list_of_c0, number_of_shifts_from_node, i_node);

                                    list_of_upper_c0 = (int*)malloc(number_of_shifts_from_node * sizeof(int));

                                    for(i_shift = 0; i_shift < number_of_shifts_from_node; i_shift++){

                                        node0 = list_of_node0[i_shift];

                                        c0 = list_of_c0[i_shift];

                                        list_of_upper_c0[i_shift] = tmp_parent_ids[node0][c0];
                                    }

                                    PRF_CSEND_INT(i_node, list_of_upper_c0, number_of_shifts_from_node, node_zero);

                                    free(list_of_node0);

                                    free(list_of_c0);

                                    free(list_of_upper_c0);
                                }
                            }

                            if(myid>0){

                                PRF_CSEND_INT(node_zero, &number_of_shifts, 1, myid);

                                list_of_node0 = (int*)malloc(number_of_shifts * sizeof(int));

                                list_of_c0 = (int*)malloc(number_of_shifts * sizeof(int));

                                //loop_shifts
                                for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                    list_of_node0[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;

                                    list_of_c0[i_shift] = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;
                                }

                                PRF_CSEND_INT(node_zero, list_of_node0, number_of_shifts, myid);

                                PRF_CSEND_INT(node_zero, list_of_c0, number_of_shifts, myid);

                                free(list_of_node0);

                                free(list_of_c0);

                                list_of_upper_c0 = (int*)malloc(number_of_shifts * sizeof(int));

                                PRF_CRECV_INT(node_zero, list_of_upper_c0, number_of_shifts, node_zero);

                                //loop_shifts
                                for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                    node0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;

                                    upper_c0 = list_of_upper_c0[i_shift];

                                    c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;

                                    upper_c1 = _C(i_layer).parent_cell[c1];

                                    w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;

                                    if(w0 == 0.0){

                                        i_tmp++;
                                    }


                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].node0 = node0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].c0 =    upper_c0;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].node1 = myid;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].c1 =    upper_c1;
                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[(i_layer + 1)][i_shift].w0 =    w0;
                                }

                                free(list_of_upper_c0);
                            }

                            if((debug_this_code)&&(i_frame == 0)){

                                report_debug(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE,
                                             "[%d] E.2: i_layer %d, number_of_shifts %d/%d w0==0 %d",
                                             myid, i_layer, number_of_shifts, number_of_shifts_in_upper_layer, i_tmp);

                            }
                        }
                    }
                }

                /* E.3 free tmp_parent_ids */
                {
                    if(tmp_parent_ids != NULL){

                        for(i_node = 0; i_node < compute_node_count; i_node++){

                            if(tmp_parent_ids[i_node] != NULL){

                                free(tmp_parent_ids[i_node]);
                            }
                        }

                        free(tmp_parent_ids);
                    }

                    if(debug_this_code){

                        report_debug(rcfd, RCFD_REPORT_CATEGORY_INITIALIZE, "[%d] E.3", myid);
                    }
                }

                i_layer++;
            }

            if(debug_this_code){

                char    filename[40];
                FILE    *fi = NULL;

                sprintf(filename,"%s_%d", "tmp_parallel_L1_shifts", myid);

                fi = fopen(filename,"w");

                i_layer = 1;

                //loop_states
                for(i_state = 0; i_state < rcfd->Solver_Dict->number_of_states; i_state++){

                    //loop_phases
                    for(i_phase = 0; i_phase < rcfd->Solver_Dict->number_of_phases; i_phase++){

                        //loop_frames
                        for(i_frame = 0; i_frame < rcfd->Solver_Dict->number_of_frames; i_frame++){

                            number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];

                            fprintf(fi, "%d\n", number_of_shifts);

                            //loop_shifts
                            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                                fprintf(fi, "%d %d %d %d %e\n", rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0, rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0,

                                    rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node1, rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1, rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0);
                            }
                        }
                    }
                }

                fclose(fi);
            }

        }
    }
#endif
#endif

#if 1 /* communication */

    void parallel_shift_data(rCFD_type *rcfd, const int i_state, const int i_layer, const int i_phase, const int i_island, const int i_frame)
    {
#if RP_NODE

#if 1   // local vars

        int     i_list, i_data, i_shift, i_node;

        int     number_of_shifts, MPI_number_of_shifts, MPI_index, c0, c1, n0, n1, index_in_total_list, loop_offset_in0, loop_offset_in1;

        double  w0, data0;

        double  *list_of_double = NULL;

#endif

        /* fill local C2Cs.shifts_out and send it to Node-0 */
        {
            if (myid > 0){

                number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_out[i_layer];

                PRF_CSEND_INT(node_zero, &number_of_shifts, 1, myid);

                if(number_of_shifts > 0){

                    MPI_number_of_shifts = number_of_shifts * rcfd->Phase_Dict[i_phase].number_of_data;

                    list_of_double = (double*)malloc(MPI_number_of_shifts * sizeof(double));

                    //loop_shifts
                    for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                        c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].c0;

                        //loop_data
                        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                            i_list = rcfd->Phase_Dict[i_phase].number_of_data * i_shift + i_data;

                            list_of_double[i_list] = _C(i_layer).data[i_phase][c0][i_data];

                            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                                w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].w0;
                                n1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].node1;

                                if((n1<0) ||(n1>node_last)){

                                    report_error(rcfd, RCFD_REPORT_CATEGORY_RUN_SHIFT, "[%d] invalid node1", myid);
                                }
                                if((w0 < 0) || (w0 > 1e5)){

                                    report_error(rcfd, RCFD_REPORT_CATEGORY_RUN_SHIFT, "[%d] invalid w0", myid);
                                }

                                rcfd->Balance[i_phase][i_data].node2node_flux[myid][n1] += w0;
                                rcfd->Balance[i_phase][i_data].node2node_data_flux[myid][n1] += w0 * _C(i_layer).data[i_phase][c0][i_data];
                            }
                        }
                    }

                    PRF_CSEND_REAL(node_zero, list_of_double, MPI_number_of_shifts, myid);

                    free(list_of_double);
                }
            }

            if (myid == 0){

                number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_out[i_layer];

                //loop_shifts
                for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                    c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].c0;

                    MPI_index = rcfd->C2Cs[i_state][i_phase][i_frame].in2out[i_layer][node_zero][i_shift];

                    //loop_data
                    for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                        rcfd->Parallel->C2Cs_MPI->shifts_in[MPI_index].data[i_data] = _C(i_layer).data[i_phase][c0][i_data];

                        if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                            w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].w0;
                            n1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_out[i_layer][i_shift].node1;

                            rcfd->Balance[i_phase][i_data].node2node_flux[myid][n1] += w0;
                            rcfd->Balance[i_phase][i_data].node2node_data_flux[myid][n1] += w0 * _C(i_layer).data[i_phase][c0][i_data];

                            if((n1<0) ||(n1>node_last)){

                                report_error(rcfd, RCFD_REPORT_CATEGORY_RUN_SHIFT, "[%d] invalid node1", myid);
                            }
                        }

                    }
                }

                /* fill rCFD_TRACER_COMM by Node>0 rCFD_TRACER_OUTGOING data */
                for(i_node = 1; i_node < compute_node_count; i_node++){

                    PRF_CRECV_INT(i_node, &number_of_shifts, 1, i_node);

                    if(number_of_shifts > 0){

                        MPI_number_of_shifts = number_of_shifts * rcfd->Phase_Dict[i_phase].number_of_data;

                        list_of_double = (double*)malloc(MPI_number_of_shifts * sizeof(double));

                        PRF_CRECV_REAL(i_node, list_of_double, MPI_number_of_shifts, i_node);

                        //loop_shifts
                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                            MPI_index = rcfd->C2Cs[i_state][i_phase][i_frame].in2out[i_layer][i_node][i_shift];

                            for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                                i_list = rcfd->Phase_Dict[i_phase].number_of_data * i_shift + i_data;

                                rcfd->Parallel->C2Cs_MPI->shifts_in[MPI_index].data[i_data] = list_of_double[i_list];
                            }
                        }

                        free(list_of_double);
                    }
                }
            }
        }

        /* Node-0 sends C2Cs.shifts_in to other Nodes; all nodes set _C(i_layer).data_shift */
        {
            if (myid == 0){

                number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_in[i_layer];

                index_in_total_list = number_of_shifts;

                for(i_node = 1; i_node < compute_node_count; i_node++){

                    number_of_shifts = rcfd->C2Cs[i_state][i_phase][i_frame].number_of_shifts_to_node_zero[i_layer][i_node];

                    PRF_CSEND_INT(i_node, &number_of_shifts, 1, node_zero);

                    if(number_of_shifts > 0){

                        MPI_number_of_shifts = number_of_shifts * rcfd->Phase_Dict[i_phase].number_of_data;

                        list_of_double = (double*)malloc(MPI_number_of_shifts * sizeof(double));

                        //loop_shifts
                        for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                            //loop_data
                            for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                                i_list = rcfd->Phase_Dict[i_phase].number_of_data * i_shift + i_data;

                                list_of_double[i_list] = rcfd->Parallel->C2Cs_MPI->shifts_in[(index_in_total_list + i_shift)].data[i_data];
                            }
                        }

                        PRF_CSEND_REAL(i_node, list_of_double, MPI_number_of_shifts, node_zero);

                        index_in_total_list += number_of_shifts;

                        free(list_of_double);
                    }


                }

                loop_offset_in0 = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[i_layer][i_island];
                loop_offset_in1 = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[i_layer][(i_island + 1)];

                for(i_shift = loop_offset_in0; i_shift < loop_offset_in1; i_shift++){

                    c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;
                    c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;
                    w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;

                    if(w0 > 0.0){

                        //loop_data
                        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                            data0 = rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].data[i_data];

                            _C(i_layer).data_shift[i_phase][c1][i_data] =

                                (w0 * data0 + _C(i_layer).weight_after_shift[c1] * _C(i_layer).data_shift[i_phase][c1][i_data])/

                                (w0 + _C(i_layer).weight_after_shift[c1]);


                            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                                n0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;

                                rcfd->Balance[i_phase][i_data].node2node_flux[n0][myid] += w0;
                                rcfd->Balance[i_phase][i_data].node2node_data_flux[n0][myid] += w0 * data0;
                            }
                        }

                        _C(i_layer).weight_after_shift[c1] += w0;
                    }
                }
            }

            if (myid > 0){

                PRF_CRECV_INT(node_zero, &number_of_shifts, 1, node_zero);

                if(number_of_shifts > 0){

                    MPI_number_of_shifts = number_of_shifts * rcfd->Phase_Dict[i_phase].number_of_data;

                    list_of_double=(double*)malloc(MPI_number_of_shifts * sizeof(double));

                    PRF_CRECV_REAL(node_zero, list_of_double, MPI_number_of_shifts, node_zero);

                    loop_offset_in0 = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[i_layer][i_island];
                    loop_offset_in1 = rcfd->C2Cs[i_state][i_phase][i_frame].island_offsets_in[i_layer][(i_island + 1)];

                    for(i_shift = loop_offset_in0; i_shift < loop_offset_in1; i_shift++){

                        c0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c0;
                        c1 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].c1;
                        w0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].w0;

                        if((w0 > 0.0)){

                            //loop_data
                            for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

                                i_list = i_shift * rcfd->Phase_Dict[i_phase].number_of_data + i_data;

                                data0 = list_of_double[i_list];

                                _C(i_layer).data_shift[i_phase][c1][i_data] =

                                    (w0 * data0 + _C(i_layer).weight_after_shift[c1] * _C(i_layer).data_shift[i_phase][c1][i_data])/

                                    (w0 + _C(i_layer).weight_after_shift[c1]);

                                if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                                    n0 = rcfd->C2Cs[i_state][i_phase][i_frame].shifts_in[i_layer][i_shift].node0;

                                    rcfd->Balance[i_phase][i_data].node2node_flux[n0][myid] += w0;

                                    rcfd->Balance[i_phase][i_data].node2node_data_flux[n0][myid] += w0 * data0;
                                }
                            }

                            _C(i_layer).weight_after_shift[c1] += w0;
                        }
                    }

                    free(list_of_double);
               }
            }
        }

#endif
    }

    void parallel_monitor_data_fluxes(rCFD_type *rcfd, const int i_state, const int i_layer, const int i_phase, const int i_island, const int i_frame)
    {
#if RP_NODE

        int     i_data;

        //loop_data
        for(i_data = 0; i_data < rcfd->Phase_Dict[i_phase].number_of_data; i_data++){

            if(rcfd->Balance_Dict[i_phase][i_data].type == RCFD_BALANCE_LOCAL){

                report_warning0(rcfd, RCFD_REPORT_CATEGORY_RUN_BALANCE,
                                "parallel_monitor_data_fluxes: Not implemented yet");
            }
        }
#endif
    }

#endif

#if 1 /* free */

    void parallel_free_all(rCFD_type *rcfd)
    {
#if RP_NODE

        int i_layer;

        if(rcfd->Parallel->MPI_Topo->MPI_Cell != NULL){

            //loop_layers
            for(i_layer = 0; i_layer < rcfd->Solver_Dict->number_of_layers; i_layer++){

                if(_MPI_Cell.number_of_ext_cells_per_node != NULL){

                    free(_MPI_Cell.number_of_ext_cells_per_node);
                }

                if(_MPI_Cell.number_of_host_cells_per_node != NULL){

                    free(_MPI_Cell.number_of_host_cells_per_node);
                }

                if(_MPI_Cell.host_of_cell != NULL){

                    free(_MPI_Cell.host_of_cell);
                }

                if(_MPI_Cell.host_of_ext_cells != NULL){

                    free(_MPI_Cell.host_of_ext_cells);
                }

                if(_MPI_Cell.hosting_cell_index != NULL){

                    free(_MPI_Cell.hosting_cell_index);
                }

                if(_MPI_Cell.host2ext_index != NULL){

                    free(_MPI_Cell.host2ext_index);
                }

                if(_MPI_Cell.data != NULL){

                    free(_MPI_Cell.data);
                }
            }

            free(rcfd->Parallel->MPI_Topo->MPI_Cell);
        }

        /*if(MPI_Topo.MPI_Face != NULL){

            //loop_layers
            for(i_layer = 0; i_layer < rcfd->Solver_Dict->number_of_layers; i_layer++){

                if(_MPI_Face.principal_face != NULL){

                    free(_MPI_Face.principal_face);
                }
            }

            free(MPI_Topo.MPI_Face);
        }*/

        free(rcfd->Parallel->MPI_Topo);

        if(rcfd->Parallel->C2Cs_MPI->shifts_in != NULL){

            int i_shift;
            int number_of_shifts = rcfd->Parallel->C2Cs_MPI->max_number_of_MPI_shifts;

            //loop_shifts
            for(i_shift = 0; i_shift < number_of_shifts; i_shift++){

                free(rcfd->Parallel->C2Cs_MPI->shifts_in[i_shift].data);
            }

            free(rcfd->Parallel->C2Cs_MPI->shifts_in);
        }

        free(rcfd->Parallel->C2Cs_MPI);
#endif
    }

#endif

/*****************************************************************************/

#if defined(RCFD_LIB_TEST_MODULES)
    int test_parallel(rCFD_type *rcfd)
    {
        int success = 0;
#if RP_NODE
        int sarray_i[3] = {0};
        int *rarray0_i = NULL;
        double sarray_r[3] = {0.};
        double *rarray0_r = NULL;

        sarray_i[0] = 3 * (myid) + 0;
        sarray_i[1] = 3 * (myid) + 1;
        sarray_i[2] = 3 * (myid) + 2;

        sarray_r[0] = 3.0 * (myid) + 0 + 0.5;
        sarray_r[1] = 3.0 * (myid) + 1 + 0.5;
        sarray_r[2] = 3.0 * (myid) + 2 + 0.5;

        PRF_GSYNC();

        Message0("\nStarting rCFD test 'parallel'\n");

        {
            int ret_val;
            int i_proc;
            int root = 0;
            int count = 2;
            int total_count = 0;
            int *srcounts   = (int*)malloc(sizeof(int) * compute_node_count);
            int *displs     = (int*)malloc(sizeof(int) * compute_node_count);

            // test integer versions of gather functions
            if(myid == root){
                rarray0_i = (int*)malloc(sizeof(int) * compute_node_count * count);
            }
            PRF_GATHER_INT(rarray0_i, sarray_i, count, root);
            if(myid == root){
                int chars = 3;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                int *expected_rarray = (int*)malloc(sizeof(int) * compute_node_count * count);
                int i_recv;
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = i_proc * count + i_recv;
                        expected_rarray[index] = 3 * (i_proc) + i_recv;
                        sprintf(&message[index*chars], "%2d ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_i[index]){
                            success |= 1;
                        }
                    }
                }
                Message("Test PRF_GATHER_INT from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = i_proc * count + i_recv;
                        sprintf(&message[index*chars], "%2d ", rarray0_i[index]);
                    }
                }
                Message("Test PRF_GATHER_INT from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_i); rarray0_i = NULL;
            PRF_GSYNC();


            total_count = 0;
            count = 3;
            for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                srcounts[i_proc] = count;
                displs[i_proc] = 3 * i_proc;
                total_count += srcounts[i_proc];
            }
            if(myid == root){
                rarray0_i = (int*)malloc(sizeof(int) * compute_node_count * 3);
            }
            PRF_GATHERV_INT(rarray0_i, sarray_i, srcounts, srcounts[myid], displs, root);
            if(myid == root){
                int chars = 3;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                int *expected_rarray = (int*)malloc(sizeof(int) * compute_node_count * count);
                int i_recv;
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = displs[i_proc] + i_recv;
                        expected_rarray[index] = 3 * (i_proc) + i_recv;
                        sprintf(&message[index*chars], "%2d ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_i[index]){
                            success |= 2;
                        }
                    }
                }
                Message("Test PRF_GATHERV_INT from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = displs[i_proc] + i_recv;
                        sprintf(&message[index*chars], "%2d ", rarray0_i[index]);
                    }
                }
                Message("Test PRF_GATHERV_INT from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_i); rarray0_i = NULL;
            PRF_GSYNC();


            srcounts[0] = srcounts[1] = 2;
            total_count = 0;
            for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                total_count += srcounts[i_proc];
            }
            ret_val = PRF_CREATE_GATHERV0_INT(&rarray0_i, sarray_i, srcounts[myid]);
            if(ret_val != total_count){
                success |= 4;
                Message("Test PRF_CREATE_GATHERV0_INT from proc %d expected %d\n", myid, total_count);
                Message("Test PRF_CREATE_GATHERV0_INT from proc %d obtained %d\n", myid, ret_val);
            }
            else if(myid == 0){
                int chars = 3;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                int *expected_rarray = (int*)malloc(sizeof(int) * compute_node_count * count);
                int i_recv;
                int index;
                for(i_proc = 0, index = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < srcounts[i_proc]; ++i_recv){
                        expected_rarray[index] = 3 * (i_proc) + i_recv;
                        sprintf(&message[index*chars], "%2d ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_i[index]){
                            success |= 4;
                        }
                        ++index;
                    }
                }
                Message("Test PRF_CREATE_GATHERV0_INT from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0, index = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < srcounts[i_proc]; ++i_recv){
                        sprintf(&message[index*chars], "%2d ", rarray0_i[index]);
                        ++index;
                    }
                }
                Message("Test PRF_CREATE_GATHERV0_INT from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_i); rarray0_i = NULL;
            PRF_GSYNC();


            // test floating point versions of gather functions
            count = 2;
            total_count = 0;

            if(myid == root){
                rarray0_r = (double*)malloc(sizeof(double) * compute_node_count * count);
            }
            PRF_GATHER_REAL(rarray0_r, sarray_r, count, root);
            if(myid == root){
                int chars = 6;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                double *expected_rarray = (double*)malloc(sizeof(double) * compute_node_count * count);
                int i_recv;
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = i_proc * count + i_recv;
                        expected_rarray[index] = 3.0 * (i_proc) + i_recv + 0.5;
                        sprintf(&message[index*chars], "%5.2f ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_r[index]){
                            success |= 1;
                        }
                    }
                }
                Message("Test PRF_GATHER_REAL from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = i_proc * count + i_recv;
                        sprintf(&message[index*chars], "%5.2f ", rarray0_r[index]);
                    }
                }
                Message("Test PRF_GATHER_REAL from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_r); rarray0_r = NULL;
            PRF_GSYNC();


            total_count = 0;
            count = 3;
            for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                srcounts[i_proc] = count;
                displs[i_proc] = 3 * i_proc;
                total_count += srcounts[i_proc];
            }
            if(myid == root){
                rarray0_r = (double*)malloc(sizeof(double) * compute_node_count * 3);
            }
            PRF_GATHERV_REAL(rarray0_r, sarray_r, srcounts, srcounts[myid], displs, root);
            if(myid == root){
                int chars = 6;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                double *expected_rarray = (double*)malloc(sizeof(double) * compute_node_count * count);
                int i_recv;
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = displs[i_proc] + i_recv;
                        expected_rarray[index] = 3.0 * (i_proc) + i_recv + 0.5;
                        sprintf(&message[index*chars], "%5.2f ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_r[index]){
                            success |= 2;
                        }
                    }
                }
                Message("Test PRF_GATHERV_REAL from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < count; ++i_recv){
                        int index = displs[i_proc] + i_recv;
                        sprintf(&message[index*chars], "%5.2f ", rarray0_r[index]);
                    }
                }
                Message("Test PRF_GATHERV_REAL from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_r); rarray0_r = NULL;
            PRF_GSYNC();


            srcounts[0] = srcounts[1] = 2;
            total_count = 0;
            for(i_proc = 0; i_proc < compute_node_count; ++i_proc){
                total_count += srcounts[i_proc];
            }
            ret_val = PRF_CREATE_GATHERV0_REAL(&rarray0_r, sarray_r, srcounts[myid]);
            if(ret_val != total_count){
                success |= 4;
                Message("Test PRF_CREATE_GATHERV0_REAL from proc %d expected %d\n", myid, total_count);
                Message("Test PRF_CREATE_GATHERV0_REAL from proc %d obtained %d\n", myid, ret_val);
            }
            else if(myid == 0){
                int chars = 6;
                char *message = (char*)malloc(sizeof(char) * compute_node_count * count * chars + 1);
                double *expected_rarray = (double*)malloc(sizeof(double) * compute_node_count * count);
                int i_recv;
                int index;
                for(i_proc = 0, index = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < srcounts[i_proc]; ++i_recv){
                        expected_rarray[index] = 3.0 * (i_proc) + i_recv + 0.5;
                        sprintf(&message[index*chars], "%5.2f ", expected_rarray[index]);
                        if(expected_rarray[index] != rarray0_r[index]){
                            success |= 4;
                        }
                        ++index;
                    }
                }
                Message("Test PRF_CREATE_GATHERV0_REAL from proc %d expected [%s]\n", myid, message);
                for(i_proc = 0, index = 0; i_proc < compute_node_count; ++i_proc){
                    for(i_recv = 0; i_recv < srcounts[i_proc]; ++i_recv){
                        sprintf(&message[index*chars], "%5.2f ", rarray0_r[index]);
                        ++index;
                    }
                }
                Message("Test PRF_CREATE_GATHERV0_REAL from proc %d received [%s]\n", myid, message);

                free(message);
                free(expected_rarray);
            }
            free(rarray0_r); rarray0_r = NULL;
            PRF_GSYNC();

            free(srcounts);
            free(displs);
        }
#endif //  RP_NODE
        return success;
    }
#endif // RCFD_LIB_TEST_MODULES

#endif // RCFD_PARALLEL
