%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                 VTOOLS                 %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%Victor Chavarrias (victor.chavarrias@deltares.nl)
%
%$Revision$
%$Date$
%$Author$
%$Id$
%$HeadURL$
%
%

function data=gdm_read_data_map(fdir_mat,fpath_map,varname,varargin)

%% PARSE

parin=inputParser;

addOptional(parin,'tim',[]);
addOptional(parin,'var_idx',[]);
addOptional(parin,'layer',[]);
addOptional(parin,'do_load',1);
addOptional(parin,'tol_t',5/60/24);
addOptional(parin,'idx_branch',[]);
addOptional(parin,'branch','');
% addOptional(parin,'bed_layers',[]); We use <layer> for flow and sediment

parse(parin,varargin{:});

time_dnum=parin.Results.tim;
var_idx=parin.Results.var_idx;
layer=parin.Results.layer;
do_load=parin.Results.do_load;
tol_t=parin.Results.tol_t;
idx_branch=parin.Results.idx_branch;
branch=parin.Results.branch;
% bed_layers=parin.Results.bed_layers;

%% CALC

% varname=D3D_var_derived2raw(varname); %I don't think I need it...
[ismor,is1d]=D3D_is(fpath_map);
[var_str,varname_changed]=D3D_var_num2str(varname,'is1d',is1d,'ismor',ismor); %in the call in <create_mat_map_2DH_01> we already change the name. Here we only need the save name. Otherwise we need to pass <simdef.D3D.structure>
% fpath_sal=mat_tmp_name(fdir_mat,var_str,'tim',time_dnum,'var_idx',var_idx,'branch',branch);
fpath_sal=mat_tmp_name(fdir_mat,var_str,'tim',time_dnum,'branch',branch); %`var_idx` cannot be in the name because it is not saved as such. 

if exist(fpath_sal,'file')==2
    if do_load
        messageOut(NaN,sprintf('Loading mat-file with raw data: %s',fpath_sal));
        load(fpath_sal,'data')
    else
        messageOut(NaN,sprintf('Mat-file with raw data exists: %s',fpath_sal));
        data=NaN;
    end
else
    messageOut(NaN,sprintf('Reading raw data for variable: %s',var_str));
    if isempty(idx_branch) %2D
        if ischar(varname)
            data=gdm_read_data_map_char(fpath_map,varname,'tim',time_dnum,'tol_t',tol_t);%,'bed_layers',bed_layers);
        else
            %outdated?
            data=gdm_read_data_map_num(fpath_map,varname,'tim',time_dnum);
        end
    else %1D
        [~,~,~,~,~,idx_tim]=D3D_time_dnum(fpath_map,time_dnum,'fdir_mat',fdir_mat);
        val=gdm_read_data_map_1D(fpath_map,var_str,idx_branch,idx_tim);
        data.val=val;
    end
    save_check(fpath_sal,'data');
end

%% layer

if ~isempty(layer)
    %maybe better to search for [layer] in the ones coming from EHY?
    idx_f=D3D_search_index_layer(data);
    data.val=submatrix(data.val,idx_f,layer);
    
end
%get desired fractions
if ~isempty(var_idx)
%     idx_f=D3D_search_index_in_dimension(data_lyrfrac,'sedimentFraction'); 
    idx_f=D3D_search_index_fraction(data); 
    data.val=submatrix(data.val,idx_f,var_idx); %take submatrix along dimension
    %sum over sediment dimension
%     data.val=sum(data.val,idx_f); %why? not here. I copied this from another location. 
end

% if ~isempty(bed_layers)
%     %maybe better to search for [layer] in the ones coming from EHY?
%    data.val=data.val(:,:,layer);
% end

end %function