 		CUDA_CHECK_RETURN( cudaMemcpy(oGRID, dev_oGRID, iMap_bytes, cudaMemcpyDeviceToHost) );
 		end_t = clock();
 		break;
+	case 3: // MEAN
+		CUDA_CHECK_RETURN( cudaMemset(dev_oGRID, 0.0,  (size_t)iMap_bytes) );
+		start_t = clock();
+		reduction_3d_mean<<<grid,block>>>( dev_iGRIDi, map_len, Nmaps, dev_oGRID );
+		CUDA_CHECK_RETURN( cudaMemcpy(oGRID, dev_oGRID, iMap_bytes, cudaMemcpyDeviceToHost) );
+		end_t = clock();
+		break;
 	}
 	// elapsed time [ms]:
 	elapsed_time += (int)( (double)(end_t  - start_t ) / (double)CLOCKS_PER_SEC * 1000 );
 
+	/*
+	 * 		****2-D aggregation****
+	 */
 	// call 2d reduction to computing mean & std for scatterplot:
 	dim3 block_2d( BLOCKSIZE*BLOCKSIZE,1,1);
 	dim3 grid_2d ( Nmaps,1,1);

		}
	}
}
__global__ void reduction_2d_sum( const double *lin_maps, unsigned int map_len, unsigned int Nmaps, double *sum_map_2d ){
	/*
	 * 		lin_maps:	|------|------|------|	...	|------|

		res					= map_len - bdx * nBlks;

		// copy the first bdx pixels to shared mem:
		sdata[tid]			= powf( lin_maps[  tid + bdx*0		+ bix*map_len ] - sum_map_2d[bix], 2 ); // deviation of pixel from mean val
		__syncthreads();

		/* each thread is responsible for its tid in every bix
		 * (excluding the residual pixels smaller than bdx):
		 */
		while( i < nBlks + remind ){
			sdata[tid]		+= powf( lin_maps[ tid + bdx*i		+ bix*map_len ] - sum_map_2d[bix], 2 ); // sum deviations
			i++;
		}
		// residual pixels not forming a block equal to bdx is handled by tid<res threads:
		if( tid<res ){
			sdata[tid]		+= powf( lin_maps[ tid + bdx*nBlks	+ bix*map_len ] - sum_map_2d[bix], 2 ); // sum deviations
		}
		__syncthreads();

