#ifndef RCFD_USER
#define RCFD_USER

#include "udf.h"
#include "rCFD_library.h"

/* (C)  2023
    Stefan Pirker
    Particulate Flow Modelling
    Johannes Kepler University, Linz, Austria
    www.particulate-flow.at
*/


    /* user Names */
    /*************************************************************************************/

#if 1

#if RP_NODE

        static short    rCFD_monitors_initiated = 0;
#endif

    // diam definitions
    
    //#define     _diam_1     (2.0 * 0.002)
    //#define     _diam_2     (2.0 * 0.001)
	
	//#define     _diam_1     (2.0 * 0.0014121)		// email ThL 12.09.23
    //#define     _diam_2     (2.0 * 0.00097609)

	#define     _diam_1     (2.0 * 0.0020375)		// email ThL 12.09.23
    #define     _diam_2     (2.0 * 0.00072931)
    
    #define     _diam_32    (2.0 * 0.001142)
	
	enum{
		_linear,
		_square,
		_square_Marco,
		_number_of_drift_approx_modes
	};
	
	#define 	_k_diam_mode _square
    
    //#define     _k_diam     	0.2
	#define     _k_diam     	0.2 				// PD2_asBD, PS4_asBD
	#define     _k_sqr_diam 	5.0e3
	#define     _k_sqrM_diam	5.0e2
    
    //#define     _m_diam_1   0.25
    //#define     _m_diam_2   0.75

	//#define     _m_diam_1   (1.0 - 0.5275)			// email ThL 12.09.23
    //#define     _m_diam_2   0.5275

	#define     _m_diam_1   (1.0 - 0.4365)			// email ThL 12.09.23
    #define     _m_diam_2   0.4365

    /* names of phases */
    enum{
        ph_gas,
        ph_solid,
        number_of_phase_names
    };

    /* names of data - phase_0 */
    enum{
        number_of_phase0_data_names
    };

    /* names of data - phase_1 */
    enum{
        c_r1,   // r1 = 0.002 (y-neg drift)
        c_r2,   // r2 = 0.001 (y-pos drift)
        
        number_of_phase1_data_names
    };

    /* names of ANSYS/fluent UDMIs */
    enum{
        c_vf_small,
        c_vf_large,
        c_vof_mix,
        c_vof,
        c_d23,
        c_vf_small_ref,            /* 5    */
        c_vf_large_ref,
        c_d23_ref,
        c_vof_average,              /* 8    */
        c_vf_small_average,
        c_vf_large_average,
        c_vf_small_average_ref,
        c_vf_large_average_ref,
        number_of_UDMI_names
    };
    
    #define     _number_of_bins     25
    
    #define     _y_min      0.00
    #define     _y_max      0.25    

#if RP_NODE

    static  int average_counter = 0;
#endif

    #define     _batch_run  0

#endif

    /* user set Dicts */
    /*************************************************************************************/

#if 1

    void user_set_solver_dict(void *handle)
    {
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "report_priority", 1);

        rcfd_set_dict_property_r(handle, "Solver_Dict", -1, -1, -1, "global_time_step", (2.0 * 5.0e-3));

        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "number_of_frames", (996 / 2));

		rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "number_of_phases", 2);

        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "face_swap_diffusion_on", 1);
		
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "balance_correction_on", 1);
		
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "data_drifting_on", 1);
		
		rcfd_set_dict_property_r(handle, "Solver_Dict", -1, -1, -1, "vof_cutoff_limit", 1.0e-6);
		
        rcfd_set_dict_property_i(handle, "Solver_Dict", -1, -1, -1, "swap_drift_set_local_mean_at_node_boundaries", 1);
    }

    void user_set_phase_dict(void *handle)
    {
#if RP_NODE

		rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_gas, -1, "number_of_data", 0);

		rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "number_of_data", 2);
		
		rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "density", 1450.0);
		
		rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "vof_max", 0.65);

		rcfd_set_dict_property_r(handle, "Phase_Dict", -1, ph_solid, -1, "d23_diam", _diam_32);
		
		rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "hindered_drift_on", 0);	

		rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "local_control_conc_sum_on", 0);	

		rcfd_set_dict_property_i(handle, "Phase_Dict", -1, ph_solid, -1, "local_drift_on", 0);
#endif
    }

    void user_set_rec_dict(void *handle)
    {
        rcfd_set_dict_property_s(handle, "Rec_Dict", -1, -1, -1, "format", "replay_format");
    }

    void user_set_data_dict(void *handle)
    {
#if RP_NODE

		double	drift_velocity[3];
		
		drift_velocity[0] = 0.0;
		
		drift_velocity[2] = 0.0;
		
		// c_r1 = large particles
		{
			rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, c_r1, "type", "concentration_data");
		
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r1, "face_swap_diffusion", 0.1);

			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r1, "d1_diam", _diam_1);
			
			rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, c_r1, "drifting_data", 1);
			
			if(_k_diam_mode == _linear) 			drift_velocity[1] = - _k_diam * (_diam_1 - _diam_32);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = - _k_sqr_diam * (_diam_1 - _diam_32) * (_diam_1 - _diam_32);	// downward drift
			else if (_k_diam_mode == _square_Marco)	drift_velocity[1] = - _k_sqrM_diam * (_diam_1*_diam_1 - _diam_32*_diam_32);			// downward drift
			else 									drift_velocity[1] = 0.0;
			
			rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r1, "drift_velocity", drift_velocity);
		}
		
		// c_r2 = small particles
		{
			rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, c_r2, "type", "concentration_data");
		
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r2, "face_swap_diffusion", 0.1);
			
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r2, "d1_diam", _diam_2);
			
			rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, c_r2, "drifting_data", 1);
			
			if(_k_diam_mode == _linear) 			drift_velocity[1] = - _k_diam * (_diam_2 - _diam_32);
			else if (_k_diam_mode == _square)		drift_velocity[1] = _k_sqr_diam * (_diam_2 - _diam_32) * (_diam_2 - _diam_32);	// upward drift
			else if (_k_diam_mode == _square_Marco)	drift_velocity[1] = - _k_sqrM_diam * (_diam_2*_diam_2 - _diam_32*_diam_32);		// upward drift
			else 									drift_velocity[1] = 0.0;
			
			rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r2, "drift_velocity", drift_velocity);
		}

#endif
    }

	void user_set_cell_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        //fluent_set_cell_dict(handle, i_layer);
		
        rcfd_set_dict_property_i(handle, "Cell_Dict", i_layer, -1, -1, "number_of_user_vars", 3);		
#endif
    }

    void user_set_face_dict(void *handle, const short i_layer)
    {
#if RP_NODE
        //fluent_set_face_dict(handle, i_layer);
#endif
    }	

#endif

    /* user init */
    /*************************************************************************************/

#if 1

    void user_init_data(void *handle, const short i_layer)
    {
#if RP_NODE

        int     i_cell;

        int     int_cell_count  = rcfd_get_int_cell_count(handle, i_layer);
        
        double  **solid_data = rcfd_get_cell_data(handle, i_layer, ph_solid);		

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

            solid_data[i_cell][c_r1] = _m_diam_1;
            
            solid_data[i_cell][c_r2] = _m_diam_2;
        }

        Message0("\n\n...user_init_Data\n");
#endif
    }

#endif

    /* user accesses Data */
    /*************************************************************************************/

#if 1

    void user_post(void *handle)
    {
#if RP_NODE

#define DEBUG_LOCAL 0

#if 1   /* local vars */

        int     i_layer, i_cell, i_bin, i_frame;
        
        int     int_cell_count;
		
		double	*number_count_small_ref = NULL, *number_count_large_ref = NULL;
        
        double	mass_of_large_part, y_mean_of_large_part, mass_of_small_part, y_mean_of_small_part, mass_of_part, y_mean_of_part;
        
        double	mass_of_large_part_ref, y_mean_of_large_part_ref, mass_of_small_part_ref, y_mean_of_small_part_ref, mass_of_part_ref, y_mean_of_part_ref;
		
		double	global_time, vol_large, vol_small, bin_width, bin_y_coord;
		
		double	*vf_small_ref = NULL, *vf_large_ref = NULL;
		
        double  *bin_CFD_large = NULL, *bin_CFD_small = NULL, *bin_rCFD_vof = NULL, *bin_rCFD_large = NULL, *bin_rCFD_small = NULL;		
		
		double  **solid_data = NULL, **vof = NULL, **x = NULL;

        Domain  *d = Get_Domain(1);

        Thread  *t, *t_mix = NULL, *t_gas;
#endif

		if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post start");
		
        // 1. access rCFD vars, init local vars
        {
            i_layer = 0;
            
			int_cell_count  = rcfd_get_int_cell_count(handle, i_layer);
			
            i_frame = rcfd_get_rec_global_frame(handle);

            solid_data = rcfd_get_cell_data(handle, i_layer, ph_solid);
			
			x = rcfd_get_cell_center(handle, i_layer); // [i_cell][i_dim]
            
            vof = rcfd_get_cell_vof(handle, i_layer, i_frame);	// [i_cell][i_phase]
			
			global_time = rcfd_get_time(handle);
			
			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post access rCFD data");
        }

        // 2. eval. reference and rCFD volume fractions
        {           
            number_count_small_ref = (double*)malloc(int_cell_count * sizeof(double));
            number_count_large_ref = (double*)malloc(int_cell_count * sizeof(double));
            
            vf_small_ref = (double*)malloc(int_cell_count * sizeof(double));
            vf_large_ref = (double*)malloc(int_cell_count * sizeof(double));

            // get particle counts from ip files, set t_mix, t_gas
            thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){

                t_mix = t;

                t_gas = THREAD_SUB_THREAD(t, 0);

                begin_c_loop_int(i_cell, t){

                    number_count_large_ref[i_cell] = C_U(i_cell, t_gas);    // ... from ip file

                    number_count_small_ref[i_cell] = C_V(i_cell, t_gas);
                    
                }end_c_loop_int(i_cell, t)
            }}
            
            // set volume fraction of small/large particles (reference)
            for(i_cell = 0; i_cell < int_cell_count; ++i_cell){

                vol_large = number_count_large_ref[i_cell] * _diam_1 * _diam_1 * _diam_1 * M_PI / 6.0;

                vol_small = number_count_small_ref[i_cell] * _diam_2 * _diam_2 * _diam_2 * M_PI / 6.0;

                if((vol_small + vol_large) > 0.0){

                    vf_small_ref[i_cell] = vol_small / (vol_small + vol_large) * vof[i_cell][ph_solid];

                    vf_large_ref[i_cell] = vol_large / (vol_small + vol_large) * vof[i_cell][ph_solid];
                }
                else{

                    vf_small_ref[i_cell] = 0.0;

                    vf_large_ref[i_cell] = 0.0;
                }
            
				// write ref volume fractions to UDMI

                C_UDMI(i_cell, t_mix, c_vf_small_ref) = vf_small_ref[i_cell];

                C_UDMI(i_cell, t_mix, c_vf_large_ref) = vf_large_ref[i_cell];
            
				// write rCFD volume fractions to UDMI

                C_UDMI(i_cell, t_mix, c_vof) = vof[i_cell][ph_solid];
                
                C_UDMI(i_cell, t_mix, c_vf_small) = solid_data[i_cell][c_r2] * vof[i_cell][ph_solid];

                C_UDMI(i_cell, t_mix, c_vf_large) = solid_data[i_cell][c_r1] * vof[i_cell][ph_solid];
            }
            
			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post eval. reference and rCFD volume fractions");
        }

        /* 3. eval y_mean_large/small (rCFD data and ref data) */
        {
            /* rCFD data */
            {

                mass_of_large_part = 0.0; y_mean_of_large_part = 0.0;
				
                mass_of_small_part = 0.0; y_mean_of_small_part = 0.0;
				
				mass_of_part = 0.0; y_mean_of_part = 0.0;
				
				for(i_cell = 0; i_cell < int_cell_count; ++i_cell){

                    mass_of_large_part += solid_data[i_cell][c_r1] * vof[i_cell][ph_solid];	// (cell_volume * solid_density) is the same for all cells ...

                    y_mean_of_large_part += x[i_cell][1] * solid_data[i_cell][c_r1] * vof[i_cell][ph_solid];

                    mass_of_small_part += solid_data[i_cell][c_r2] * vof[i_cell][ph_solid];

                    y_mean_of_small_part += x[i_cell][1] * solid_data[i_cell][c_r2] * vof[i_cell][ph_solid];

                    mass_of_part += vof[i_cell][ph_solid];

                    y_mean_of_part += x[i_cell][1] * vof[i_cell][ph_solid];
                }

                mass_of_large_part =    PRF_GRSUM1(mass_of_large_part);

                y_mean_of_large_part =  PRF_GRSUM1(y_mean_of_large_part);

                mass_of_small_part =    PRF_GRSUM1(mass_of_small_part);

                y_mean_of_small_part =  PRF_GRSUM1(y_mean_of_small_part);

                mass_of_part =          PRF_GRSUM1(mass_of_part);

                y_mean_of_part =        PRF_GRSUM1(y_mean_of_part);

                if(mass_of_large_part > 0.0){

                    y_mean_of_large_part /= mass_of_large_part;
                }

                if(mass_of_small_part > 0.0){

                    y_mean_of_small_part /= mass_of_small_part;
                }

                if(mass_of_part > 0.0){

                    y_mean_of_part /= mass_of_part;
                }
            }

            /* ref data */
            {

                mass_of_large_part_ref = 0.0; y_mean_of_large_part_ref = 0.0;
				
                mass_of_small_part_ref = 0.0; y_mean_of_small_part_ref = 0.0;
				
				mass_of_part_ref = 0.0; y_mean_of_part_ref = 0.0;
				
				for(i_cell = 0; i_cell < int_cell_count; ++i_cell){

                    mass_of_large_part_ref += number_count_large_ref[i_cell] * _diam_1 * _diam_1 * _diam_1 * M_PI / 6.0;

                    y_mean_of_large_part_ref += x[i_cell][1] * number_count_large_ref[i_cell] * _diam_1 * _diam_1 * _diam_1 * M_PI / 6.0;


                    mass_of_small_part_ref += number_count_small_ref[i_cell] * _diam_2 * _diam_2 * _diam_2 * M_PI / 6.0;

                    y_mean_of_small_part_ref += x[i_cell][1] * number_count_small_ref[i_cell] * _diam_2 * _diam_2 * _diam_2 * M_PI / 6.0;


                    mass_of_part_ref += M_PI / 6.0 * 

                        (number_count_large_ref[i_cell] * _diam_1 * _diam_1 * _diam_1 + number_count_small_ref[i_cell] * _diam_2 * _diam_2 * _diam_2);

                    y_mean_of_part_ref += x[i_cell][1] *  M_PI / 6.0 *

                        (number_count_large_ref[i_cell] * _diam_1 * _diam_1 * _diam_1 + number_count_small_ref[i_cell] * _diam_2 * _diam_2 * _diam_2);
                }

                mass_of_large_part_ref =    PRF_GRSUM1(mass_of_large_part_ref);

                y_mean_of_large_part_ref =  PRF_GRSUM1(y_mean_of_large_part_ref);

                mass_of_small_part_ref =    PRF_GRSUM1(mass_of_small_part_ref);

                y_mean_of_small_part_ref =  PRF_GRSUM1(y_mean_of_small_part_ref);

                mass_of_part_ref =          PRF_GRSUM1(mass_of_part_ref);

                y_mean_of_part_ref =        PRF_GRSUM1(y_mean_of_part_ref);

                if(mass_of_large_part_ref > 0.0){

                    y_mean_of_large_part_ref /= mass_of_large_part_ref;
                }

                if(mass_of_small_part_ref > 0.0){

                    y_mean_of_small_part_ref /= mass_of_small_part_ref;
                }

                if(mass_of_part_ref > 0.0){

                    y_mean_of_part_ref /= mass_of_part_ref;
                }
            }
			
			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post eval. y_mean");
        }

        /* 4. node-0 writes y-mean values to monitor file */
        if(myid == 0){

            FILE    *f_out = NULL;

            if(rCFD_monitors_initiated == 0){

                if(_k_diam_mode == _linear) f_out = fopen("./monitor_rCFD.out","w");
				if(_k_diam_mode == _square) f_out = fopen("./monitor_sqr_rCFD.out","w");
				if(_k_diam_mode == _square_Marco) f_out = fopen("./monitor_sqrM_rCFD.out","w");

                rCFD_monitors_initiated = 1;
            }
            else{
                if(_k_diam_mode == _linear) f_out = fopen("./monitor_rCFD.out","a");
				if(_k_diam_mode == _square) f_out = fopen("./monitor_sqr_rCFD.out","a");
				if(_k_diam_mode == _square_Marco) f_out = fopen("./monitor_sqrM_rCFD.out","a");
            }

            if(f_out == NULL){

                rcfd_report0(handle, 0, 2, "ERROR: Could not open monitor file\n");

                return;
            }

            fprintf(f_out, "%e %e %e %e %e %e %e %e %e %e %e\n", global_time,

                y_mean_of_part, y_mean_of_large_part, y_mean_of_small_part,

                y_mean_of_part_ref, y_mean_of_large_part_ref, y_mean_of_small_part_ref,

                mass_of_part_ref, mass_of_large_part_ref, mass_of_small_part_ref, mass_of_part);

            fclose(f_out);
        }

		if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post writes y-mean values to monitor file");
		
        /* 5. average in between 4s and 5s (based on UDMI values) */
        {
            if((global_time > 4.0) && (global_time < 5.0)){

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

                    // rCFD
                    
                    C_UDMI(i_cell, t_mix, c_vof_average) = (double)average_counter * C_UDMI(i_cell, t_mix, c_vof_average) +
                    
                        C_UDMI(i_cell, t_mix, c_vof);

                    C_UDMI(i_cell, t_mix, c_vof_average) /= (double)(average_counter + 1);
                    
                    C_UDMI(i_cell, t_mix, c_vf_small_average) = (double)average_counter * C_UDMI(i_cell, t_mix, c_vf_small_average) + 
                    
                        C_UDMI(i_cell, t_mix, c_vf_small);

                    C_UDMI(i_cell, t_mix, c_vf_small_average) /= (double)(average_counter + 1);
                    
                    C_UDMI(i_cell, t_mix, c_vf_large_average) = (double)average_counter * C_UDMI(i_cell, t_mix, c_vf_large_average) + 
                    
                        C_UDMI(i_cell, t_mix, c_vf_large);

                    C_UDMI(i_cell, t_mix, c_vf_large_average) /= (double)(average_counter + 1);

                    // CFD
                    
                    C_UDMI(i_cell, t_mix, c_vf_small_average_ref) = (double)average_counter * C_UDMI(i_cell, t_mix, c_vf_small_average_ref) + 
                    
                        C_UDMI(i_cell, t_mix, c_vf_small_ref);
                    
                    C_UDMI(i_cell, t_mix, c_vf_small_average_ref) /= (double)(average_counter + 1);                 

                    C_UDMI(i_cell, t_mix, c_vf_large_average_ref) = (double)average_counter * C_UDMI(i_cell, t_mix, c_vf_large_average_ref) + 
                    
                        C_UDMI(i_cell, t_mix, c_vf_large_ref);

                    C_UDMI(i_cell, t_mix, c_vf_large_average_ref) /= (double)(average_counter + 1);
                }

                average_counter++;
            }
            else{
                
                C_UDMI(i_cell, t_mix, c_vof_average) = 0.0;
                
                C_UDMI(i_cell, t_mix, c_vf_small_average) = 0.0;
                C_UDMI(i_cell, t_mix, c_vf_large_average) = 0.0;
                
                C_UDMI(i_cell, t_mix, c_vf_small_average_ref) = 0.0;
                C_UDMI(i_cell, t_mix, c_vf_large_average_ref) = 0.0;    
            }
			
			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post average");
        }

        /* 6. write vertical profile of average conc */
        {
            // generate vertical bins
            {
                bin_rCFD_vof = (double*)malloc(_number_of_bins * sizeof(double));

                bin_rCFD_small = (double*)malloc(_number_of_bins * sizeof(double));

                bin_rCFD_large = (double*)malloc(_number_of_bins * sizeof(double));
                
                bin_CFD_small = (double*)malloc(_number_of_bins * sizeof(double));
                
                bin_CFD_large = (double*)malloc(_number_of_bins * sizeof(double));

                for(i_bin = 0; i_bin < _number_of_bins; i_bin ++){
                    
                    bin_rCFD_vof[i_bin] = 0.0;
                    
                    bin_rCFD_small[i_bin] = 0.0;    
                    
                    bin_rCFD_large[i_bin] = 0.0;
                    
                    bin_CFD_small[i_bin] = 0.0;
                    
                    bin_CFD_large[i_bin] = 0.0;     
                }
            }
                
            // loop cells and fill bins (CFD, rCFD)
            {
                for(i_cell = 0; i_cell < int_cell_count; ++i_cell){
                    
                    if((x[i_cell][0] > 0.09) && (x[i_cell][0] < 0.11)){
                    
                        i_bin = (int) floor((x[i_cell][1] - _y_min) / (_y_max - _y_min) * (double) _number_of_bins);
                        
                        if((i_bin >= 0) && (i_bin < _number_of_bins)){
                            
                            bin_rCFD_vof[i_bin] += C_UDMI(i_cell, t_mix, c_vof_average);
                            
                            bin_rCFD_small[i_bin] += C_UDMI(i_cell, t_mix, c_vf_small_average) * C_UDMI(i_cell, t_mix, c_vof_average);
                            
                            bin_rCFD_large[i_bin] += C_UDMI(i_cell, t_mix, c_vf_large_average) * C_UDMI(i_cell, t_mix, c_vof_average);
                            
                            bin_CFD_small[i_bin] += C_UDMI(i_cell, t_mix, c_vf_small_average_ref) * C_UDMI(i_cell, t_mix, c_vof_average);
                            
                            bin_CFD_large[i_bin] += C_UDMI(i_cell, t_mix, c_vf_large_average_ref) * C_UDMI(i_cell, t_mix, c_vof_average);
                        }
                    }
                }
            }
            
            // global reduction
            {
                for(i_bin = 0; i_bin < _number_of_bins; i_bin ++){
                    
                    bin_rCFD_vof[i_bin] = PRF_GRSUM1(bin_rCFD_vof[i_bin]);
                    
                    bin_rCFD_small[i_bin] = PRF_GRSUM1(bin_rCFD_small[i_bin]);  
                    
                    bin_rCFD_large[i_bin] = PRF_GRSUM1(bin_rCFD_large[i_bin]);
                    
                    bin_CFD_small[i_bin] = PRF_GRSUM1(bin_CFD_small[i_bin]);
                    
                    bin_CFD_large[i_bin] = PRF_GRSUM1(bin_CFD_large[i_bin]);    
                }               
            }

            // mean vf's
            {
                
                for(i_bin = 0; i_bin < _number_of_bins; i_bin ++){
                    
                    if(bin_rCFD_vof[i_bin] > 0.0){
                    
                        bin_rCFD_small[i_bin] /= bin_rCFD_vof[i_bin];   
                        
                        bin_rCFD_large[i_bin] /= bin_rCFD_vof[i_bin];
                    
                        bin_CFD_small[i_bin]  /= bin_rCFD_vof[i_bin];   
                        
                        bin_CFD_large[i_bin]  /= bin_rCFD_vof[i_bin];
                    }
                }               
            }
            
            // node-0 writes profiles in "w" mode
            if(myid == 0){

                FILE    *f_out = NULL;

                if(_k_diam_mode == _linear) f_out = fopen("./profiles.out","w");
				if(_k_diam_mode == _square) f_out = fopen("./profiles_sqr.out","w");
				if(_k_diam_mode == _square_Marco) f_out = fopen("./profiles_sqrM.out","w");

                if(f_out == NULL){

                    Message0("\nERROR: Could not open profile file\n");

                    return;
                }
                
                for(i_bin = 0; i_bin < _number_of_bins; i_bin ++){
                    
                    bin_width = (double)( _y_max - _y_min) / (double)_number_of_bins;
                    
                    bin_y_coord = _y_min + i_bin * bin_width + bin_width / 2.0;
                    
                    fprintf(f_out, "%e %e %e %e %e \n", 
                    
                        bin_y_coord, bin_rCFD_large[i_bin], bin_rCFD_small[i_bin], bin_CFD_large[i_bin], bin_CFD_small[i_bin]);
                }

                fclose(f_out);
            }

			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post write vertical profiles");
        
        }
        
        /* 7. free local vars */
        {
            free(number_count_large_ref);
            free(number_count_small_ref);

            free(vf_large_ref);
            free(vf_small_ref);
            
            free(bin_rCFD_vof);
            free(bin_rCFD_large);
            free(bin_rCFD_small);
            
            free(bin_CFD_large);
            free(bin_CFD_small);

			if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post free");
        }

        if(DEBUG_LOCAL) rcfd_report0(handle, 0, 2, "user_post end");

#undef DEBUG_LOCAL
		
#endif
    }

#endif

    /* register user functions */
    /*************************************************************************************/

    void user_register_functions(void *handle)
    {
        REGISTER_RCFD_UDF( user_set_solver_dict );
        REGISTER_RCFD_UDF( user_set_phase_dict );
        REGISTER_RCFD_UDF( user_set_rec_dict );
        REGISTER_RCFD_UDF( user_set_data_dict );

        REGISTER_RCFD_UDF( user_init_data );
        REGISTER_RCFD_UDF( user_post );
		
        REGISTER_RCFD_UDF( user_set_cell_dict );
        REGISTER_RCFD_UDF( user_set_face_dict );
    }

#endif
