function [fh2osfc,qrunoff,qdrai,fsat,fh2osfcall,qrunoffall,qdraiall,zwt,zwtall,twsall,et,etall,rain] = ...
         read_elm(datadir, tag, ntot, numc, ind2)
    
    strs = strsplit(tag,'_');
    ex = strs{2};
    
    if strcmp(strs{1},'SUS')
        files = dir([datadir tag '/Susquehanna_half*elm.h0.*.nc']);
    else
        files = dir([datadir tag '/Pantanal_half*elm.h0.*.nc']);
    end
    fh2osfc    = NaN(ntot,length(files));
    qrunoff    = NaN(ntot,length(files));
    qdrai      = NaN(ntot,length(files));
    fsat       = NaN(numc,length(files));
    fh2osfcall = NaN(numc*ntot,length(files));
    qrunoffall = NaN(numc*ntot,length(files));
    qdraiall = NaN(numc*ntot,length(files));
    zwtall   = NaN(numc*ntot,length(files));
    twsall   = NaN(numc*ntot,length(files));
    et    = NaN(ntot,length(files));
    etall = NaN(numc*ntot,length(files));
    rain  = NaN(numc,length(files));
    
    for i = 1 : length(files)
        filename = fullfile(files(i).folder,files(i).name);
        disp(['i = ' num2str(i) ': ' filename]);
        data1 = ncread(filename,'FH2OSFC');
        rainall = ncread(filename,'RAIN');
        rain(:,i) = rainall(1:numc);
        fh2osfcall(:,i) = data1;
        for j = 1 : ntot
            tmp = data1((j-1)*numc+1:j*numc,:);
            fh2osfc(j,i) = nanmean(tmp,1);
    %         if j == ind1(1)
    %             fh2osfcall(:,i) = tmp;
    %         end
        end

        data2 = ncread(filename,'EFLX_LH_TOT');
        data3 = ncread(filename,'QRUNOFF');
        data4 = ncread(filename,'FSAT');
        data5 = ncread(filename,'QDRAI');
        data6 = ncread(filename,'ZWT');
        %data7 = ncread(filename,'TWS');
        fsat(:,i) = data4(1:numc);
        etall(:,i)      = data2;
        qrunoffall(:,i) = data3;
        qdraiall(:,i)   = data5;
        zwtall(:,i)    = data6;
        %twsall(:,i)    = data7;
        for j = 1 : ntot
            tmp = data2((j-1)*numc+1:j*numc,:);
            tmp2= data3((j-1)*numc+1:j*numc,:);
            tmp5= data5((j-1)*numc+1:j*numc,:);
            tmp6= data6((j-1)*numc+1:j*numc,:);
            et(j,i) = nanmean(tmp,1);
            qrunoff(j,i) = nanmean(tmp2,1);
            qdrai(j,i)   = nanmean(tmp5,1);
            zwt(j,i)     = nanmean(tmp6,1);
        end
    end

end

