%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of NS2DDV.                                                      %
%                                                                                   %
% Copyright(C) 2011-2018    C. Calgaro  (caterina.calgaro@math.univ-lille1.fr)      %
%                           E. Creusé   (emmanuel.creuse@math.univ-lille1.fr)       %
%                           T. Goudon   (thierry.goudon@inria.fr)                   %
%                           A. Mouton   (alexandre.mouton@math.univ-lille1.fr)      %
%                                                                                   %
% NS2DDV is free software: you can redistribute it and/or modify it under the terms %
% of the GNU General Public License as published by the Free Software Foundation,   %
% either version 3 of the License, or (at your option) any later version.           %
%                                                                                   %
% NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY         %
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   %
% PARTICULAR PURPOSE. See the GNU General Public License for more details.          %
%                                                                                   %
% You should have received a copy of the GNU General Public License along with      %
% NS2DDV. If not, see <http://www.gnu.org/licenses/>.                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[] = plot_insitu(outputs, MESH, time, HANDLES, PARAMETERS)

	% INPUT
	% outputs			Structure gathering all the required data for visualization and saving
	% MESH				Mesh metadata
	% time				Time value
	% HANDLES			Metadata for in-situ visualization
	% PARAMETERS		List of simulation parameters

	if max(strcmp(PARAMETERS.OUTPUT.XDISPLAY, {'SINGLE_FIGURE', 'MULTIPLE_FIGURE'}))
		XDISPLAY_CASES = '''MESH'', ''SUBTRI_U''';
		XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''P'', ''U_VECTOR'', ''UX'', ''UY'', ');
		XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''VORTIC'', ''STREAM'', ''U_MODULE'', ');
		XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''UX_DX'', ''UX_DY'', ''UY_DX'', ''UY_DY'', ''SHEAR''');
		if PARAMETERS.FE.BC_NATURAL
			XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''W_VECTOR'', ''WX'', ''WY'', ''R''');
		end
		if strcmp(PARAMETERS.MODEL, 'NSDV')
			XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''RHO''');
		end
		if PARAMETERS.TESTCASE_ANALYTIC
			XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''P_EX'', ''U_VECTOR_EX'', ''UX_EX'', ''UY_EX'', ');
			XDISPLAY_CASES = strcat(XDISPLAY_CASES, '''VORTIC_EX'', ''STREAM_EX'', ''U_MODULE_EX'', ');
			XDISPLAY_CASES = strcat(XDISPLAY_CASES, ...
						'''UX_DX_EX'', ''UX_DY_EX'', ''UY_DX_EX'', ''UY_DY_EX'', ''SHEAR_EX''');
			if strcmp(PARAMETERS.MODEL, 'NSDV')
				XDISPLAY_CASES = strcat(XDISPLAY_CASES, ', ''RHO_EX''');
			end
		end

		if strcmp(PARAMETERS.OUTPUT.XDISPLAY, 'SINGLE_FIGURE')
			figure(HANDLES.h);
			globtitle = sprintf('Time t = %e', time);
			meshtitle = HANDLES.meshtitle;
			subtriutitle = sprintf('%s %s', PARAMETERS.FE.TYPE, HANDLES.subtriutitle);
			uvecttitle = HANDLES.uvecttitle;
			uxtitle = HANDLES.uxtitle;
			uytitle = HANDLES.uytitle;
			dxuxtitle = HANDLES.dxuxtitle;
			dyuxtitle = HANDLES.dyuxtitle;
			dxuytitle = HANDLES.dxuytitle;
			dyuytitle = HANDLES.dyuytitle;
			sheartitle = HANDLES.sheartitle;
			ptitle = HANDLES.ptitle;
			dxptitle = HANDLES.dxptitle;
			dyptitle = HANDLES.dyptitle;
			gradptitle = HANDLES.gradptitle;
			omegatitle = HANDLES.omegatitle;
			psititle = HANDLES.psititle;
			umodtitle = HANDLES.umodtitle;
			rhoextitle = HANDLES.rhoextitle;
			if PARAMETERS.FE.BC_NATURAL
				wvecttitle = HANDLES.wvecttitle;
				wxtitle = HANDLES.wxtitle;
				wytitle = HANDLES.wytitle;
				rtitle = HANDLES.rtitle;
			end
			if strcmp(PARAMETERS.MODEL, 'NSDV')
				rhotitle = HANDLES.rhotitle;
			end
			if PARAMETERS.TESTCASE_ANALYTIC
				uvectextitle = HANDLES.uvectextitle;
				uxextitle = HANDLES.uxextitle;
				uyextitle = HANDLES.uyextitle;
				dxuxextitle = HANDLES.dxuxextitle;
				dyuxextitle = HANDLES.dyuxextitle;
				dxuyextitle = HANDLES.dxuyextitle;
				dyuyextitle = HANDLES.dyuyextitle;
				shearextitle = HANDLES.shearextitle;
				pextitle = HANDLES.pextitle;
				dxpextitle = HANDLES.dxptitle;
				dypextitle = HANDLES.dyptitle;
				gradpextitle = HANDLES.gradptitle;
				omegaextitle = HANDLES.omegaextitle;
				psiextitle = HANDLES.psiextitle;
				umodextitle = HANDLES.umodextitle;
				if strcmp(PARAMETERS.MODEL, 'NSDV')
					rhoextitle = HANDLES.rhoextitle;
				end
			end

		elseif strcmp(PARAMETERS.OUTPUT.XDISPLAY, 'MULTIPLE_FIGURE')
			meshtitle = sprintf('%s at time = %e', HANDLES.meshtitle, A.time);
			subtriutitle = sprintf('%s %s at time = %e', A.METHOD_U, HANDLES.subtriutitle, A.time);
			uvecttitle = sprintf('%s at time = %e', HANDLES.uvecttitle, A.time);
			uxtitle = sprintf('%s at time = %e', HANDLES.uxtitle, A.time);
			uytitle = sprintf('%s at time = %e', HANDLES.uytitle, A.time);
			dxuxtitle = sprintf('%s at time = %e', HANDLES.dxuxtitle, A.time);
			dyuxtitle = sprintf('%s at time = %e', HANDLES.dyuxtitle, A.time);
			dxuytitle = sprintf('%s at time = %e', HANDLES.dxuytitle, A.time);
			dyuytitle = sprintf('%s at time = %e', HANDLES.dyuytitle, A.time);
			sheartitle = sprintf('%s at time = %e', HANDLES.sheartitle, A.time);
			ptitle = sprintf('%s at time = %e', HANDLES.ptitle, A.time);
			dxptitle = sprintf('%s at time = %e', HANDLES.dxptitle, A.time);
			dyptitle = sprintf('%s at time = %e', HANDLES.dyptitle, A.time);
			gradptitle = sprintf('%s at time = %e', HANDLES.gradptitle, A.time);
			omegatitle = sprintf('%s at time = %e', HANDLES.omegatitle, A.time);
			psititle = sprintf('%s at time = %e', HANDLES.psititle, A.time);
			umodtitle = sprintf('%s at time = %e', HANDLES.umodtitle, A.time);

			if PARAMETERS.FE.BC_NATURAL
				rtitle = sprintf('%s at time = %e', HANDLES.rtitle, A.time);
				wvecttitle = sprintf('%s at time = %e', HANDLES.wvecttitle, A.time);
				wxtitle = sprintf('%s at time = %e', HANDLES.wxtitle, A.time);
				wytitle = sprintf('%s at time = %e', HANDLES.wytitle, A.time);
			end
			if strcmp(PARAMETERS.MODEL, 'NSDV')
				rhotitle = sprintf('%s at time = %e', HANDLES.rhotitle, A.time);
			end
			if PARAMETERS.TESTCASE_ANALYTIC
				uvectextitle = sprintf('%s at time = %e', HANDLES.uvectextitle, A.time);
				uxextitle = sprintf('%s at time = %e', HANDLES.uxextitle, A.time);
				uyextitle = sprintf('%s at time = %e', HANDLES.uyextitle, A.time);
				dxuxextitle = sprintf('%s at time = %e', HANDLES.dxuxextitle, A.time);
				dyuxextitle = sprintf('%s at time = %e', HANDLES.dyuxextitle, A.time);
				dxuyextitle = sprintf('%s at time = %e', HANDLES.dxuyextitle, A.time);
				dyuyextitle = sprintf('%s at time = %e', HANDLES.dyuyextitle, A.time);
				shearextitle = sprintf('%s at time = %e', HANDLES.shearextitle, A.time);
				pextitle = sprintf('%s at time = %e', HANDLES.pextitle, A.time);
				dxpextitle = sprintf('%s at time = %e', HANDLES.dxpextitle, A.time);
				dypextitle = sprintf('%s at time = %e', HANDLES.dypextitle, A.time);
				gradpextitle = sprintf('%s at time = %e', HANDLES.gradpextitle, A.time);
				omegaextitle = sprintf('%s at time = %e', HANDLES.omegaextitle, A.time);
				psiextitle = sprintf('%s at time = %e', HANDLES.psiextitle, A.time);
				umodextitle = sprintf('%s at time = %e', HANDLES.umodextitle, A.time);
				if strcmp(PARAMETERS.MODEL, 'NSDV')
					rhoextitle = sprintf('%s at time = %e', HANDLES.rhoextitle, A.time);
				end
			end
		end
	end

	switch PARAMETERS.MODEL
		case {'NS'}
			[levels_p, levels_u, levels_umod, levels_omega, levels_psi] = isovalues_ns(time, PARAMETERS);
			levels_rho = [];
		case {'NSDV'}
			[levels_rho, levels_p, levels_u, levels_umod, levels_omega, levels_psi] = ...
				isovalues_nsdv(time, PARAMETERS);
		otherwise
			error('Unknown model');
	end
	levels_dxp = [];
	levels_dyp = [];
	levels_ux = [];
	levels_uy = [];
	levels_dxux = [];
	levels_dyux = [];
	levels_dxuy = [];
	levels_dyuy = [];
	levels_shear = [];
	levels_wx = [];
	levels_wy = [];
	levels_r = [];


	% Prepare the isovalues for scalar diagnostics
	ind_rho = [];
	ind_p = [];
	ind_dxp = [];
	ind_dyp = [];
	ind_ux = [];
	ind_uy = [];
	ind_omega = [];
	ind_psi = [];
	ind_shear = [];
	ind_umod = [];
	ind_dxux = [];
	ind_dyux = [];
	ind_dxuy = [];
	ind_dyuy = [];
	ind_wx = [];
	ind_wy = [];
	ind_r = [];
	for i=1:size(outputs,1)
		if max(strcmp(outputs{i,2}, {'PLOT', 'PLOTSAVE', 'SAVEPLOT'}))
			% We must plot the data

			switch outputs{i,1}
				case {'MESH', 'SUBTRI_U', 'U_VECTOR', 'U_VECTOR_EX', 'W_VECTOR', 'GRAD_P'}
					% Nothing to do

				case {'RHO', 'RHO_EX'}
					ind_rho = [ind_rho, i];
					levels_rho = compute_levels_insitu(outputs, ind_rho, levels_rho, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'P', 'P_EX'}
					ind_p = [ind_p, i];
					levels_p = compute_levels_insitu(outputs, ind_p, levels_p, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'P_DX', 'P_DX_EX'}
					ind_dxp = [ind_dxp, i];
					levels_dxp = compute_levels_insitu(outputs, ind_dxp, levels_dxp, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'P_DY', 'P_DY_EX'}
					ind_dyp = [ind_dyp, i];
					levels_dyp = compute_levels_insitu(outputs, ind_dyp, levels_dyp, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UX', 'UX_EX'}
					ind_ux = [ind_ux, i];
					levels_ux = compute_levels_insitu(outputs, ind_ux, levels_ux, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UY', 'UY_EX'}
					ind_uy = [ind_uy, i];
					levels_uy = compute_levels_insitu(outputs, ind_uy, levels_uy, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'U_MODULE', 'U_MODULE_EX'}
					ind_umod = [ind_umod, i];
					levels_umod = compute_levels_insitu(outputs, ind_umod, levels_umod, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'VORTIC', 'VORTIC_EX'}
					ind_omega = [ind_omega, i];
					levels_omega = compute_levels_insitu(outputs, ind_omega, levels_omega, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'STREAM', 'STREAM_EX'}
					ind_psi = [ind_psi, i];
					levels_psi = compute_levels_insitu(outputs, ind_psi, levels_psi, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UX_DX', 'UX_DX_EX'}
					ind_dxux = [ind_dxux, i];
					levels_dxux = compute_levels_insitu(outputs, ind_dxux, levels_dxux, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UX_DY', 'UX_DY_EX'}
					ind_dyux = [ind_dyux, i];
					levels_dyux = compute_levels_insitu(outputs, ind_dyux, levels_dyux, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UY_DX', 'UY_DX_EX'}
					ind_dxuy = [ind_dxuy, i];
					levels_dxuy = compute_levels_insitu(outputs, ind_dxuy, levels_dxuy, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'UY_DY', 'UY_DY_EX'}
					ind_dyuy = [ind_dyuy, i];
					levels_dyuy = compute_levels_insitu(outputs, ind_dyuy, levels_dyuy, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'SHEAR', 'SHEAR_EX'}
					ind_shear = [ind_shear, i];
					levels_shear = compute_levels_insitu(outputs, ind_shear, levels_shear, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'WX'}
					ind_wx = [ind_wx, i];
					levels_wx = compute_levels_insitu(outputs, ind_wx, levels_wx, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'WY'}
					ind_wy = [ind_wy, i];
					levels_wy = compute_levels_insitu(outputs, ind_wy, levels_wy, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				case {'R'}
					ind_r = [ind_r, i];
					levels_r = compute_levels_insitu(outputs, ind_r, levels_r, ...
									PARAMETERS.TESTCASE_ANALYTIC);

				otherwise
					error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
							outputs{i,1}, XDISPLAY_CASES));
			end
		end
	end

	for i=1:size(outputs,1)
		if max(strcmp(outputs{i,2}, {'PLOT', 'PLOTSAVE', 'SAVEPLOT'}))
			% We must plot the data

			% Identify handle, levels, title, data
			switch outputs{i,1}
				case {'MESH'}
					h = HANDLES.hmesh;
					ti = meshtitle;

				case {'SUBTRI_U'}
					h = HANDLES.hsubtriu;
					ti = subtriutitle;

				case {'RHO'}
					h = HANDLES.hrho;
					levels = levels_rho;
					ti = rhotitle;
					data = outputs{i,3};

				case {'RHO_EX'}
					h = HANDLES.hrhoex;
					levels = levels_rho;
					ti = rhoextitle;
					data = outputs{i,3};

				case {'P'}
					h = HANDLES.hp;
					levels = levels_p;
					ti = ptitle;
					data = outputs{i,3};

				case {'P_EX'}
					h = HANDLES.hpex;
					levels = levels_p;
					ti = pextitle;
					data = outputs{i,3};

				case {'P_DX'}
					h = HANDLES.hdxp;
					levels = levels_dxp;
					ti = dxptitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,1);

				case {'P_DX_EX'}
					h = HANDLES.hdxpex;
					levels = levels_dxpex;
					ti = dxpextitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,1);

				case {'P_DY'}
					h = HANDLES.hdyp;
					levels = levels_dyp;
					ti = dyptitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,2);

				case {'P_DY_EX'}
					h = HANDLES.hdypex;
					levels = levels_dypex;
					ti = dypextitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,2);

				case {'GRAD_P'}
					h = HANDLES.hgradp;
					ti = gradptitle;
					data = outputs{i,3};

				case {'GRAD_P_EX'}
					h = HANDLES.hgradpex;
					ti = gradpextitle;
					data = outputs{i,3};

				case {'U_VECTOR'}
					h = HANDLES.huvect;
					ti = uvecttitle;
					data = outputs{i,3};

				case {'U_VECTOR_EX'}
					h = HANDLES.huvectex;
					ti = uvectextitle;
					data = outputs{i,3};

				case {'UX'}
					h = HANDLES.hux;
					levels = levels_ux;
					ti = uxtitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,1);

				case {'UX_EX'}
					h = HANDLES.huxex;
					levels = levels_ux;
					ti = uxextitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,1);

				case {'UY'}
					h = HANDLES.huy;
					levels = levels_uy;
					ti = uytitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,2);

				case {'UY_EX'}
					h = HANDLES.huy_ex;
					levels = levels_uy;
					ti = uyextitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,2);

				case {'UX_DX'}
					h = HANDLES.hdxux;
					levels = levels_dxux;
					ti = dxuxtitle;
					data = outputs{i,3};

				case {'UX_DX_EX'}
					h = HANDLES.hdxuxex;
					levels = levels_dxux;
					ti = dxuxextitle;
					data = outputs{i,3};

				case {'UX_DY'}
					h = HANDLES.hdyux;
					levels = levels_dyux;
					ti = dyuxtitle;
					data = outputs{i,3};

				case {'UX_DY_EX'}
					h = HANDLES.hdyuxex;
					levels = levels_dyux;
					ti = dyuxextitle;
					data = outputs{i,3};

				case {'UY_DX'}
					h = HANDLES.hdxuy;
					levels = levels_dxuy;
					ti = dxuytitle;
					data = outputs{i,3};

				case {'UY_DX_EX'}
					h = HANDLES.hdxuyex;
					levels = levels_dxuy;
					ti = dxuyextitle;
					data = outputs{i,3};

				case {'UY_DY'}
					h = HANDLES.hdyuy;
					levels = levels_dyuy;
					ti = dyuytitle;
					data = outputs{i,3};

				case {'UY_DY_EX'}
					h = HANDLES.hdyuyex;
					levels = levels_dyuy;
					ti = dyuyextitle;
					data = outputs{i,3};

				case {'SHEAR'}
					h = HANDLES.hshear;
					levels = levels_shear;
					ti = sheartitle;
					data = outputs{i,3};

				case {'SHEAR_EX'}
					h = HANDLES.hshearex;
					levels = levels_shear;
					ti = shearextitle;
					data = outputs{i,3};

				case {'VORTIC'}
					h = HANDLES.homega;
					levels = levels_omega;
					ti = omegatitle;
					data = outputs{i,3};

				case {'VORTIC_EX'}
					h = HANDLES.homegaex;
					levels = levels_omega;
					ti = omegaextitle;
					data = outputs{i,3};

				case {'STREAM'}
					h = HANDLES.hpsi;
					levels = levels_psi;
					ti = psititle;
					data = outputs{i,3};

				case {'STREAM_EX'}
					h = HANDLES.hpsiex;
					levels = levels_psi;
					ti = psiextitle;
					data = outputs{i,3};

				case {'U_MODULE'}
					h = HANDLES.humod;
					levels = levels_umod;
					ti = umodtitle;
					data = outputs{i,3};

				case {'U_MODULE_EX'}
					h = HANDLES.humodex;
					levels = levels_umod;
					ti = umodextitle;
					data = outputs{i,3};

				case {'W_VECTOR'}
					h = HANDLES.hwvect;
					ti = wvecttitle;
					data = outputs{i,3};

				case {'WX'}
					h = HANDLES.hwx;
					levels = levels_wx;
					ti = wxtitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,1);

				case {'WY'}
					h = HANDLES.hwy;
					levels = levels_uy;
					ti = wytitle;
					% Warning: we must redirect data to the corresponding vector field
					ivect = outputs{i,3}(1);
					data = outputs{ivect,3}(:,2);

				case {'R'}
					h = HANDLES.hr;
					levels = levels_r;
					ti = rtitle;
					data = outputs{i,3};

				otherwise
					error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
							outputs{i,1}, XDISPLAY_CASES));
			end

			% Identify the mesh
			switch outputs{i,1}
				case {'MESH', 'R', 'RHO', 'RHO_EX', 'P', 'P_EX', 'P_DX', 'P_DX_EX', 'P_DY', 'P_DY_EX', ...
							'VORTIC', 'VORTIC_EX', 'STREAM', 'STREAM_EX', ...
							'UX_DX', 'UX_DX_EX', 'UX_DY', 'UX_DY_EX', 'UY_DX', 'UY_DX_EX', ...
							'UY_DY', 'UY_DY_EX', 'SHEAR', 'SHEAR_EX'}
					p = MESH.p1;
					t = MESH.t1(1:3,:);

				case {'GRAD_P', 'GRAD_P_EX'}
					p = zeros(2,size(MESH.t1,2));
					for ke=1:size(p,2)
						p(:,ke) = (MESH.p1(:,MESH.t1(1,ke)) + MESH.p1(:,MESH.t1(2,ke)) + ...
										MESH.p1(:,MESH.t1(3,ke)))/3.;
					end

				case {'SUBTRI_U', 'U_VECTOR', 'U_VECTOR_EX', 'UX', 'UX_EX', 'UY', 'UY_EX', ...
							'U_MODULE', 'U_MODULE_EX', 'W_VECTOR', 'WX', 'WY'}
					switch PARAMETERS.FE.TYPE
						case {'P1B'}
							p = MESH.p1b;
							t = MESH.st1b;
						case {'P2'}
							p = MESH.p2;
							t = MESH.st2;
						otherwise
							error('Unknown FE discretization');
					end

				otherwise
					error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
							outputs{i,1}, XDISPLAY_CASES));
			end

			% Plot data
			switch_plot(h, PARAMETERS.OUTPUT.XDISPLAY);
			switch outputs{i,1}
				case {'MESH', 'SUBTRI_U'}
					% Plot a mesh
					plot_mesh(p, t, ti, h);

				case {'R', 'RHO', 'RHO_EX', 'P', 'P_EX', 'VORTIC', 'VORTIC_EX', 'STREAM', 'STREAM_EX', ...
							'UX_DX', 'UX_DY', 'UY_DX', 'UY_DY', 'SHEAR', ...
							'UX_DX_EX', 'UX_DY_EX', 'UY_DX_EX', 'UY_DY_EX', 'SHEAR_EX', ...
							'UX', 'UY', 'U_MODULE', 'UX_EX', 'UY_EX', 'U_MODULE_EX', 'WX', 'WY'}
					plot_contour(data, p, t, levels, ti, h);
					%plot_color(data, p, t, ti, h);

				case {'P_DX', 'P_DX_EX', 'P_DY', 'P_DY_EX'}
					plot_color(data, p, t, ti, h);

				case {'U_VECTOR', 'U_VECTOR_EX', 'W_VECTOR', 'GRAD_P', 'GRAD_P_EX'}
					plot_vector(data, p, ti, h);

				otherwise
					error(sprintf('Cannot plot the diagnostic %s\nOnly the following diagnostics are allowed: %s', ...
							outputs{i,1}, XDISPLAY_CASES));
			end
			if strcmp(PARAMETERS.OUTPUT.XDISPLAY, 'SINGLE_FIGURE')
				suplabel(globtitle, 't');
				drawnow;
			end
		end
	end


	%%%%%%%%%%%%%%%%%%%%%%%
	% Auxiliary functions %
	%%%%%%%%%%%%%%%%%%%%%%%
	function[newlevels] = compute_levels_insitu(outputs, ind, oldlevels, TEST_ANALYTIC)
		newlevels = [];
		if TEST_ANALYTIC
			% We wait for the reading of both data and exact data
			if (numel(ind) == 2)
				if (numel(oldlevels) > 0)
					% We already have the isovalues --> nothing to do
					newlevels = oldlevels;
				else
					% No pre-built isovalues --> compute them
					isomin = min(min([outputs{ind(1),3},outputs{ind(1),3}]));
					isomax = max(max([outputs{ind(1),3},outputs{ind(1),3}]));
					newlevels = linspace(isomin,isomax,7);
				end
			end
		else
			if (numel(oldlevels) > 0)
				% We already have the isovalues --> nothing to do
				newlevels = oldlevels;
			else
				% No pre-built isovalues --> compute them
				isomin = min(min(outputs{ind(1),3}));
				isomax = max(max(outputs{ind(1),3}));
				newlevels = linspace(isomin,isomax,7);
			end
		end
	end

end