function [ntot, q, fdfall, fdf, R2, NSE, KGE, ind1, ind2, imax] = read_mosart(datadir, tag, yr1, yr2,mu_gsim, lon_gsim,lat_gsim,yr, area, ntot2)
    
    strs = strsplit(tag,'_');
    ex = strs{2};
    
    if strcmp(ex,'all') || strcmp(ex,'run3') || strcmp(ex,'run4') || strcmp(ex,'run5') 
        ntot = 1000;
    elseif strcmp(tag,'outputs_geo_CLMCRUNCEPv7') || strcmp(tag,'outputs_geo_CLMGSWP3v1') || strcmp(tag,'outputs_geo_noinund_CLMGSWP3v1')
        ntot = 100;
    elseif strcmp(tag,'outputs_sub_CLMCRUNCEPv7') || strcmp(tag,'outputs_sub_CLMGSWP3v1')
        ntot = 100;
    elseif strcmp(ex,'run2') || strcmp(ex,'sur3') || strcmp(ex,'sur2')
        ntot = 100;
    elseif strcmp(ex,'run') || strcmp(ex,'sur') 
        ntot = 500;
    elseif strcmp(ex,'default') 
        if strcmp(strs{1},'SUS')
            ntot = 1;
        else
            ntot = 10;
        end
    elseif strcmp(ex,'man') || strcmp(ex,'man2')
        ntot = 50;
    elseif strcmp(ex,'hydrology')
        ntot = 200;
    elseif strcmp(ex,'hydraulic') 
        ntot = 200;
    elseif strcmp(ex,'hydraulic2')
        ntot = 1000;
    else
        ntot = ntot2;
    end
    
    if strcmp(strs{1},'SUS')
        files = dir([datadir tag '/Susquehanna_half_*.mosart.h0.*.nc']);
        if length(files) == 0
            files = dir([datadir tag '/Pantanal_half_*.mosart.h0.*.nc']);
        end
    else
        files = dir([datadir tag '/Pantanal_half_*.mosart.h0.*.nc']);
    end
    q  = NaN(ntot,length(files));
    %q2 = NaN(ntot,length(files));
    fdf = NaN(ntot,length(files));
    
    for i = 1 : length(files)
        filename = fullfile(files(i).folder,files(i).name);

        disp(files(i).name);

        q1(:,1) = ncread(filename,'RIVER_DISCHARGE_OVER_LAND_LIQ');
        q1(:,2) = ncread(filename,'RIVER_DISCHARGE_TO_OCEAN_LIQ');
        if i == 1
            numc = size(q1,1)/ntot;
            fdfall = NaN(numc*ntot,length(files));
        end
        q1 = nanmean(q1,2);
        if i == 1
            if strcmp(strs{1},'SUS')
            [ioutlet,icontribution]   = find_mosart_cell('inputdata_run/MOSART_Susquehanna_run3_template_c230518.nc',lon_gsim,lat_gsim,area);    
            else
            [ioutlet,icontribution]   = find_mosart_cell('inputdata_run/MOSART_Pantanal_all_template_c220126.nc',lon_gsim,lat_gsim,area);
            end
            %[ioutlet2,icontribution2] = find_mosart_cell('inputdata_run/MOSART_Pantanal_all_template_c220126.nc',lon_gsim2,lat_gsim2,area2);
        end
        q(:,i)  = q1(ioutlet : numc : end);
        %q2(:,i) = q1(ioutlet2: numc : end);

        if ~contains(tag,'noinund') && ~contains(ex,'default')
            fdfall(:,i) = ncread(filename,'FLOODED_FRACTION');
        end
    end
    
    for i = 1 : ntot
        [R2(i),~,NSE(i)] = estimate_evaluation_metric(mu_gsim(yr >= yr1 & yr <= yr2),q(i,:)');
        KGE(i) = estimateKGE(mu_gsim(yr >= yr1 & yr <= yr2),q(i,:)');
    end
    ind1 = find(NSE == max(NSE));
    ind2 = find(R2  == max(R2));

    k = 1;
    imax = NaN(ntot,1);
    for i = 1 : ntot
        tmp = nanmean(reshape(q(i,:)',[12,yr2-yr1+1]),2);
        [~,imax(i)] = max(tmp);
        k = k + 1;
    end
end