#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

    #define     _diam_23    (2.0 * 0.001142)
    
	enum{
		_linear,
		_square,
		_number_of_drift_approx_modes
	};

	#define 	_k_diam_mode _square

    #define     _k_diam     	0.2
	#define     _k_sqr_diam 	5.0e3

    #define     _number_of_particle_fractions   5

    #define     _diam_1     (2. * _diam_23)
    #define     _diam_2     (3./2. * _diam_23)
    #define     _diam_3     _diam_23
    #define     _diam_4     (3./4. * _diam_23)
    #define     _diam_5     (1./2. * _diam_23)

    #define     _mf_diam_1  0.15
    #define     _mf_diam_2  0.2
    #define     _mf_diam_3  0.427
    #define     _mf_diam_4  0.123
    #define     _mf_diam_5  0.1

    #define     loop_diam   for(i_diam = 0; i_diam < _number_of_particle_fractions; i_diam ++)
    
    static short    rCFD_monitors_initiated =   0;
	static short	rCFD_profiles_initiated =   0;
	static short    rCFD_y_mean_bed_initiated = 0;
    
#endif

    /* 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.002284 (y-neg drift)
        c_r2,   // r2 = 0.001713
        c_r3,   // r3 = 0.001142 (y-pos drift)
        c_r4,   // r4 = 0.0008565 
        c_r5,   // r5 = 0.000571 (y-pos drift)
        
        _number_of_diams
    };

    enum{
        c_vf_diam_1,
        c_vf_diam_2,
        c_vf_diam_3,
        c_vf_diam_4,
        c_vf_diam_5,

        c_vf_diam_1_ref,
        c_vf_diam_2_ref,
        c_vf_diam_3_ref,
        c_vf_diam_4_ref,
        c_vf_diam_5_ref,

        c_vof_ref,				//10
        
        c_vf_diam_1_av,
        c_vf_diam_2_av,
        c_vf_diam_3_av,
        c_vf_diam_4_av,
        c_vf_diam_5_av,

        c_vf_diam_1_ref_av,
        c_vf_diam_2_ref_av,
        c_vf_diam_3_ref_av,
        c_vf_diam_4_ref_av,
        c_vf_diam_5_ref_av,

        number_of_UDMI_names
    };

    #define     _number_of_bins     25
    
    #define     _y_min      0.00
    #define     _y_max      0.25
    
    #define     _y_bottom   0.06
	
	#define		_y_bottom_line	0.005	// 0.005

    
#if RP_NODE
    
    static  double  part_diam[5];   // set by rCFD_user_init_Data

    static  int     average_counter = 0;

#endif

	#define		loop_bins	for(i_bin = 0; i_bin < _number_of_bins; i_bin++)
    
#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", (300 / 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", 5);
		
		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_23);
		
		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", "quarter_jumps_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, "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_23);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = - _k_sqr_diam * (_diam_1 - _diam_23) * (_diam_1 - _diam_23);	// 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 = large 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, "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_23);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = - _k_sqr_diam * (_diam_2 - _diam_23) * (_diam_2 - _diam_23);	// downward 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);
		}

		// c_r3 = d23 particles
		{
			rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, c_r3, "type", "concentration_data");
			
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r3, "d1_diam", _diam_3);
			
			rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, c_r3, "drifting_data", 1);
			
			if(_k_diam_mode == _linear) 			drift_velocity[1] = - _k_diam * (_diam_3 - _diam_23);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = - _k_sqr_diam * (_diam_3 - _diam_23) * (_diam_3 - _diam_23);	// no drift
			else 									drift_velocity[1] = 0.0;
			
			rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r3, "drift_velocity", drift_velocity);
		}

		// c_r4 = small particles
		{
			rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, c_r4, "type", "concentration_data");
			
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r4, "d1_diam", _diam_4);
			
			rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, c_r4, "drifting_data", 1);
			
			if(_k_diam_mode == _linear) 			drift_velocity[1] = - _k_diam * (_diam_4 - _diam_23);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = _k_sqr_diam * (_diam_4 - _diam_23) * (_diam_4 - _diam_23);	// upward drift
			else 									drift_velocity[1] = 0.0;
			
			rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r4, "drift_velocity", drift_velocity);
		}

		// c_r5 = Small particles
		{
			rcfd_set_dict_property_s(handle, "Data_Dict", -1, ph_solid, c_r5, "type", "concentration_data");
			
			rcfd_set_dict_property_r(handle, "Data_Dict", -1, ph_solid, c_r5, "d1_diam", _diam_5);
			
			rcfd_set_dict_property_i(handle, "Data_Dict", -1, ph_solid, c_r5, "drifting_data", 1);
			
			if(_k_diam_mode == _linear) 			drift_velocity[1] = - _k_diam * (_diam_5 - _diam_23);	
			else if (_k_diam_mode == _square)		drift_velocity[1] = _k_sqr_diam * (_diam_5 - _diam_23) * (_diam_5 - _diam_23);	// upward drift
			else 									drift_velocity[1] = 0.0;
			
			rcfd_set_dict_property_r_3d(handle, "Data_Dict", -1, ph_solid, c_r5, "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] = _mf_diam_1;            
            solid_data[i_cell][c_r2] = _mf_diam_2;
            solid_data[i_cell][c_r3] = _mf_diam_3;            
            solid_data[i_cell][c_r4] = _mf_diam_4;
            solid_data[i_cell][c_r5] = _mf_diam_5;            
        }

        part_diam[0] = _diam_1;
        part_diam[1] = _diam_2;
        part_diam[2] = _diam_3;
        part_diam[3] = _diam_4;
        part_diam[4] = _diam_5;
		
        Message0("\n\n...user_init_Data\n");
#endif
    }

#endif

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

#if 1

    void user_post(void *handle)
    {
#if RP_NODE

#if 1   /* local vars */
        int     i_layer, i_cell, i_data, i_frame, i_diam, i_bin;
        
        int     int_cell_count; 
        
        double  vol_part_ref_all, vf_all_ref, vol_part, vof_max;
        
        double  vol_bottom, vol_top, vol_bottom_ref, vol_top_ref, global_time, sum_of_denom, d32_bottom;
		
		double  part_vol, acc_part_vol, acc_part_vol_y, y_mean_bed;
        
        double  *vol_bottom_diam = NULL, *vol_top_diam = NULL, *vol_bottom_diam_ref = NULL, *vol_top_diam_ref = NULL;
        
        double  *mf_bottom_diam = NULL, *mf_top_diam = NULL, *mf_bottom_diam_ref = NULL, *mf_top_diam_ref = NULL;
		
		double	*cell_volume = NULL, **cell_center = NULL;
		   
        double  **number_count_ref = NULL, number_count_ref_all;
		
		double	**x = NULL, **solid_data = NULL, **vof = NULL;
		
		double	*bin_rCFD_vol = NULL, *bin_ref_vol = NULL, **bin_rCFD_diam = NULL, **bin_ref_diam = NULL, bin_width, bin_y_coord;

        Domain  *d = Get_Domain(1);

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

        // access vars, set 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]
			
			vof_max = 0.65;
			
			cell_volume = rcfd_get_cell_volume(handle, i_layer);
			
			cell_center = rcfd_get_cell_center(handle, i_layer);

			global_time = rcfd_get_time(handle);
		}
        
        // get number_count_ref[i_diam]
        {
            number_count_ref = (double**)malloc(_number_of_particle_fractions * sizeof(double*));
            
            loop_diam{
                
                number_count_ref[i_diam] = (double*)malloc(int_cell_count * sizeof(double));
                
                for(i_cell = 0; i_cell < int_cell_count; ++i_cell){
                    
                    number_count_ref[i_diam][i_cell] = 0.0;
                }
            }

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

                t_solid = THREAD_SUB_THREAD(t, 1);

                begin_c_loop_int(i_cell, t){

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

                    number_count_ref[1][i_cell] = C_V(i_cell, t_gas);

                    number_count_ref[2][i_cell] = C_W(i_cell, t_gas);

                    number_count_ref[3][i_cell] = C_U(i_cell, t_solid);

                    number_count_ref[4][i_cell] = C_V(i_cell, t_solid);
                    
                }end_c_loop_int(i_cell, t)
            }}          
        }
        
        // set UDMI(vf_diam_ref)
        {
            thread_loop_c(t,d){if(FLUID_CELL_THREAD_P(t)){begin_c_loop_int(i_cell, t){

                number_count_ref_all = 0.0;
                
                vol_part_ref_all = 0.0;
                
                loop_diam{
                    
                    number_count_ref_all += number_count_ref[i_diam][i_cell];
                    
                    vol_part_ref_all += number_count_ref[i_diam][i_cell] * 
                    
                        part_diam[i_diam] * part_diam[i_diam] * part_diam[i_diam] * M_PI / 6.0;
                }
                
                if(vol_part_ref_all < (cell_volume[i_cell] * vof_max)){
                    
                    vf_all_ref = vol_part_ref_all / cell_volume[i_cell];
                }
                else{
                    
                    vf_all_ref = vof_max;
                }           
                
                C_UDMI(i_cell, t, c_vof_ref) = vf_all_ref;
                
                if(number_count_ref_all > 0.0){
                    
                    loop_diam{
                        
                        vol_part = number_count_ref[i_diam][i_cell] * 
                        
                            part_diam[i_diam] * part_diam[i_diam] * part_diam[i_diam] * M_PI / 6.0;
                        
                        C_UDMI(i_cell, t, (c_vf_diam_1_ref + i_diam)) = vol_part / vol_part_ref_all * vf_all_ref;
                    }
                }
                else{
                    
                    loop_diam{
                        
                        C_UDMI(i_cell, t, (c_vf_diam_1_ref + i_diam)) = 0.0;
                    }
                }

            }end_c_loop_int(i_cell, t)}}                        
        }
        
        // set UDMI(vf_diam)
        {
            for(i_cell = 0; i_cell < int_cell_count; ++i_cell){
                
                for(i_data = 0; i_data < _number_of_particle_fractions; i_data++){
                    
                    C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_data)) = solid_data[i_cell][i_data] * vof[i_cell][ph_solid];
                }
            }
        }
        
        // set UDMI(vf_diam_av), UDMI(vf_diam_ref_av)
        {
            if((global_time > 4.0) && (global_time < 5.0)){

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

                    loop_diam{
                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_av + i_diam)) = 
                        
                            (double)average_counter * C_UDMI(i_cell, t_mix, (c_vf_diam_1_av + i_diam)) + 
                            
                            C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_diam));
                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_av + i_diam)) /= (double)(average_counter + 1);

                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref_av + i_diam)) = 
                        
                            (double)average_counter * C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref_av + i_diam)) + 
                            
                            C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref + i_diam));
                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref_av + i_diam)) /= (double)(average_counter + 1);
                    }
                }

                average_counter++;
            }
            else{
                
                for(i_cell = 0; i_cell < int_cell_count; ++i_cell){

                    loop_diam{
                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_av + i_diam)) = 0.0;
                        
                        C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref_av + i_diam)) = 0.0;
                    }
                }
            }
            
        }
        
		// set and write y_mean_bed
		{
			acc_part_vol_y = 0.0;
			acc_part_vol =   0.0;
			
			for(i_cell = 0; i_cell < int_cell_count; i_cell++){
				
				for(i_diam = 0; i_diam < _number_of_diams; i_diam++){
					
					part_vol = part_diam[i_diam] * part_diam[i_diam] * part_diam[i_diam] * M_PI / 6.0;
					
				    acc_part_vol_y += number_count_ref[i_diam][i_cell] * part_vol * cell_center[i_cell][1];
					
					acc_part_vol += number_count_ref[i_diam][i_cell] * part_vol;
				}
			}

			acc_part_vol_y = PRF_GRSUM1(acc_part_vol_y);
			acc_part_vol   = PRF_GRSUM1(acc_part_vol);

			if(acc_part_vol > 0.0) y_mean_bed = acc_part_vol_y / acc_part_vol;
			else 				   y_mean_bed = 0.0;

            // node-0 writes y_mean_bed.out
            if(myid == 0){

                FILE    *f_out = NULL;

                if(rCFD_y_mean_bed_initiated == 0){

					f_out = fopen("./y_mean_bed.out","w");

                    rCFD_y_mean_bed_initiated = 1;
                }
                else

					f_out = fopen("./y_mean_bed.out","a");

                fprintf(f_out, "%e %e\n", global_time, y_mean_bed);

                fclose(f_out);
            }

		}
		
        // write vertical profiles
        {
            // generate vertical bins
            {
                bin_rCFD_vol = (double*)malloc(_number_of_bins * sizeof(double));
                bin_ref_vol = (double*)malloc(_number_of_bins * sizeof(double));

                bin_rCFD_diam = (double**)malloc(_number_of_bins * sizeof(double*));
                bin_ref_diam = (double**)malloc(_number_of_bins * sizeof(double*));
				
				loop_bins{
					
					bin_rCFD_vol[i_bin] = 0.0;
					bin_ref_vol[i_bin] = 0.0;
					
					bin_rCFD_diam[i_bin] = (double*)malloc(_number_of_particle_fractions * sizeof(double));
					bin_ref_diam[i_bin] = (double*)malloc(_number_of_particle_fractions * sizeof(double));
					
					loop_diam{
						
						bin_rCFD_diam[i_bin][i_diam] = 0.0;
						bin_ref_diam[i_bin][i_diam] =  0.0;
					}
				}
            }
                
            // loop cells and fill bins (ref, 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_vol[i_bin] += vof[i_cell][ph_solid] * cell_volume[i_cell];
							
							bin_ref_vol[i_bin] += C_UDMI(i_cell, t_mix, c_vof_ref) * cell_volume[i_cell];
                            
                            loop_diam{
								
								bin_rCFD_diam[i_bin][i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_diam)) * vof[i_cell][ph_solid] * cell_volume[i_cell];
 								
								bin_ref_diam[i_bin][i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref + i_diam)) * C_UDMI(i_cell, t_mix, c_vof_ref) * cell_volume[i_cell];
							}
                        }
                    }
                }
            }
            
            // global reduction
            {
                loop_bins{
                    
                    bin_rCFD_vol[i_bin] = PRF_GRSUM1(bin_rCFD_vol[i_bin]);
                    bin_ref_vol[i_bin] = PRF_GRSUM1(bin_ref_vol[i_bin]);
                    
                    loop_diam{
						
						bin_rCFD_diam[i_bin][i_diam] = PRF_GRSUM1(bin_rCFD_diam[i_bin][i_diam]);  
						bin_ref_diam[i_bin][i_diam] = PRF_GRSUM1(bin_ref_diam[i_bin][i_diam]); 
                    }
                }               
            }

            // mean vf's
            {   
                loop_bins{
                    
                    if(bin_rCFD_vol[i_bin] > 0.0){
                    
                        loop_diam{
							
							bin_rCFD_diam[i_bin][i_diam] /= bin_rCFD_vol[i_bin];   
						}
                    }

                    if(bin_ref_vol[i_bin] > 0.0){
                    
                        loop_diam{
							
							bin_ref_diam[i_bin][i_diam] /= bin_ref_vol[i_bin];   
						}
					}
                }               
            }
            
            // node-0 writes profiles
            if(myid == 0){

                FILE    *f_out = NULL;
				
                if(rCFD_profiles_initiated == 0){

					if(_k_diam_mode == _linear) f_out = fopen("./profiles.out","w");
				    if(_k_diam_mode == _square) f_out = fopen("./profiles_sqr.out","w");
					
                    rCFD_profiles_initiated = 1;
                }
                else{
					if(_k_diam_mode == _linear) f_out = fopen("./profiles.out","a");
				    if(_k_diam_mode == _square) f_out = fopen("./profiles_sqr.out","a");
				}

                if(f_out == NULL){

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

                    return;
                }
                
                loop_bins{
                    
                    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 ", bin_y_coord);

					loop_diam{
						
						fprintf(f_out, "%e ", bin_rCFD_diam[i_bin][i_diam]);
					}

					loop_diam{
						
						fprintf(f_out, "%e ", bin_ref_diam[i_bin][i_diam]);
					}
					
					fprintf(f_out, "\n");
                }

                fclose(f_out);
            }
        }
        
        // write bottom/top psd
        {
            // bottom / top masses (rCFD, ref)
            {
                // initializations
                {
                    vol_bottom      = 0.0;              
                    vol_top         = 0.0;
                    
                    vol_bottom_ref  = 0.0;              
                    vol_top_ref     = 0.0;
                    
                    vol_bottom_diam =       (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    vol_bottom_diam_ref =   (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    vol_top_diam =          (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    vol_top_diam_ref =      (double*)malloc( _number_of_particle_fractions * sizeof(double));
                    
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] =       0.0;
                        vol_bottom_diam_ref[i_diam] =   0.0;
                        vol_top_diam[i_diam] =          0.0;
                        vol_top_diam_ref[i_diam] =      0.0;
                    }
                }
                
                // accumulated fractional bottom/top volumes
				for(i_cell = 0; i_cell < int_cell_count; ++i_cell){
                    
                    if(x[i_cell][1] < _y_bottom){
                        
                        // rCFD
                        vol_bottom += vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_bottom_diam[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_diam)) * cell_volume[i_cell];
                        }
                        
                        // ref
                        vol_bottom_ref += C_UDMI(i_cell, t_mix, c_vof_ref) * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_bottom_diam_ref[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref + i_diam)) * cell_volume[i_cell];
                        }
                    }
                    else{
                        
                        // rCFD
                        vol_top += vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_top_diam[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_diam)) * cell_volume[i_cell];
                        }
                        
                        // ref
                        vol_top_ref += C_UDMI(i_cell, t_mix, c_vof_ref) * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_top_diam_ref[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1_ref + i_diam)) * cell_volume[i_cell];
                        }
                    }
                }
            
                // global reduction
                {
                    vol_bottom =        PRF_GRSUM1(vol_bottom);
                    vol_bottom_ref =    PRF_GRSUM1(vol_bottom_ref);
                    vol_top =           PRF_GRSUM1(vol_top);
                    vol_top_ref =       PRF_GRSUM1(vol_top_ref);
                    
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] =       PRF_GRSUM1(vol_bottom_diam[i_diam]);
                        vol_bottom_diam_ref[i_diam] =   PRF_GRSUM1(vol_bottom_diam_ref[i_diam]);
                        vol_top_diam[i_diam] =          PRF_GRSUM1(vol_top_diam[i_diam]);
                        vol_top_diam_ref[i_diam] =      PRF_GRSUM1(vol_top_diam_ref[i_diam]);
                    }   
                }
                
                // mf_diam top/bottom
                {
                    // init mf's
                    {
                        mf_bottom_diam =        (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        mf_top_diam =           (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        mf_bottom_diam_ref =    (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        mf_top_diam_ref =       (double*)malloc( _number_of_particle_fractions * sizeof(double));
                        
                        loop_diam{
                            
                            mf_bottom_diam[i_diam] =        0.0;
                            mf_top_diam[i_diam] =           0.0;
                            mf_bottom_diam_ref[i_diam] =    0.0;
                            mf_top_diam_ref[i_diam] =       0.0;
                        }
                    }
                    
                    if(vol_bottom > 0.0){
                        
                        loop_diam{
                            
                            mf_bottom_diam[i_diam] = vol_bottom_diam[i_diam] / vol_bottom;
                        }
                    }
                    
                    if(vol_top > 0.0){
                        
                        loop_diam{
                            
                            mf_top_diam[i_diam] = vol_top_diam[i_diam] / vol_top;
                        }
                    }
                    
                    if(vol_bottom_ref > 0.0){
                        
                        loop_diam{
                            
                            mf_bottom_diam_ref[i_diam] = vol_bottom_diam_ref[i_diam] / vol_bottom_ref;
                        }
                    }
                    
                    if(vol_top_ref > 0.0){
                        
                        loop_diam{
                            
                            mf_top_diam_ref[i_diam] = vol_top_diam_ref[i_diam] / vol_top_ref;
                        }
                    }
                }
            
			}
            
            // write monitor file
            if(myid == 0){

                FILE    *f_out = NULL;

                if(rCFD_monitors_initiated == 0){

					if(_k_diam_mode == _linear) f_out = fopen("./rCFD_monitor_poly.out","w");
					if(_k_diam_mode == _square) f_out = fopen("./rCFD_monitor_sqr_poly.out","w");

                    rCFD_monitors_initiated = 1;
                }
                else{
					
					if(_k_diam_mode == _linear) f_out = fopen("./rCFD_monitor_poly.out","a");
					if(_k_diam_mode == _square) f_out = fopen("./rCFD_monitor_sqr_poly.out","a");
                }

                if(f_out == NULL){

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

                    return;
                }

                fprintf(f_out, "%e ", global_time);
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_bottom_diam[i_diam]);
                }
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_top_diam[i_diam]);
                }
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_bottom_diam_ref[i_diam]);
                }
                
                loop_diam{
                    
                    fprintf(f_out, "%e ", mf_top_diam_ref[i_diam]);
                }
                
                fprintf(f_out, "\n");

                fclose(f_out);
            }
        
        }
        
        // write bottom line d32
        {
            // bottom masses (rCFD, ref)
            {
                // initializations
                {
                    vol_bottom      = 0.0;              
                    
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] =       0.0;
                    }
                }
                
                // accumulated fractional bottom/top volumes
				for(i_cell = 0; i_cell < int_cell_count; ++i_cell){
                    
                    if(x[i_cell][1] < _y_bottom_line){
                        
                        // rCFD
                        vol_bottom += vof[i_cell][ph_solid] * cell_volume[i_cell];
                        
                        loop_diam{
                            
                            vol_bottom_diam[i_diam] += C_UDMI(i_cell, t_mix, (c_vf_diam_1 + i_diam)) * cell_volume[i_cell];
                        }
                    }
                }
            
                // global reduction
                {
                    vol_bottom =        PRF_GRSUM1(vol_bottom);
                    
                    loop_diam{
                        
                        vol_bottom_diam[i_diam] = PRF_GRSUM1(vol_bottom_diam[i_diam]);
                    }   
                }
                
                // mf_diam bottom
                {
                    // init mf's
					{
						
						loop_diam{
							
							mf_bottom_diam[i_diam] = 0.0;
						}
					}
                     
                    if(vol_bottom > 0.0){
                        
                        loop_diam{
                            
                            mf_bottom_diam[i_diam] = vol_bottom_diam[i_diam] / vol_bottom;
                        }
                    }
                }
				
				// d32_bottom
				{
					
					sum_of_denom = 0.0;
					
					loop_diam{
						
						sum_of_denom += mf_bottom_diam[i_diam] / part_diam[i_diam];
					}
					
					if (sum_of_denom > 0.0){
						
						d32_bottom = 1.0 / sum_of_denom;
					}
					else{
						
						d32_bottom = 0.0;
					}
				}
			
			}
            
            // write monitor file
            if(myid == 0){

                FILE    *f_out = NULL;

                if(rCFD_monitors_initiated == 0){

					if(_k_diam_mode == _linear) f_out = fopen("./rCFD_monitor_poly_d32.out", "w");
					if(_k_diam_mode == _square) f_out = fopen("./rCFD_monitor_sqr_poly_d32.out", "w");

                    rCFD_monitors_initiated = 1;
                }
                else{
					if(_k_diam_mode == _linear) f_out = fopen("./rCFD_monitor_poly_d32.out", "a");
					if(_k_diam_mode == _square) f_out = fopen("./rCFD_monitor_sqr_poly_d32.out", "a");
                }

                if(f_out == NULL){

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

                    return;
                }

                fprintf(f_out, "%e %e\n", global_time, d32_bottom);

                fclose(f_out);
            }
        
        }
    
        // free local vars
        {
            loop_diam{
                
                free(number_count_ref[i_diam]);
            }
            
            free(number_count_ref);
    
            free(vol_bottom_diam);
            free(vol_bottom_diam_ref);
            free(vol_top_diam);
            free(vol_top_diam_ref);
    
            free(mf_bottom_diam);
            free(mf_bottom_diam_ref);
            free(mf_top_diam);
            free(mf_top_diam_ref);
			
			loop_bins{
				
				free(bin_rCFD_diam[i_bin]);
				free(bin_ref_diam[i_bin]);
			}
			
			free(bin_rCFD_diam);
			free(bin_ref_diam);
        }
        
        Message0("\n\n...user_post\n");
#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
