function [XXi, FFdi,YY, CtyCons, gammas, VAcoeff, alphas, Surplusn, VAn, Pi, PiF,XXi_gs,FFdi_gs] = WIOTcleaning_f0515_DC(FFread, XXread, VVread, tauX, tauF, J, S, N)
    %0515第一次修改。增加关税额的计算。
    %cleaning Data
    VV = VVread;
    VV(VV<0) = 0;%后面就再也没用过了，注意是否还会再用

    FF = FFread;
    NN = zeros(J,1); 
    for i=1:N
        NN(:,1) = NN(:,1) + FF(:,(i*5));%表示这一行所有存货相加的和
        FF1(:,i) = sum(FF(:,(((i-1)*5)+1):(i*5)),2);
        
    end
    XX = XXread;
    
    YY = sum(XX,2)+sum(FF1,2);%sum后为2意思是计算行向量的值。
   
    %check1 = YY-Ya; 
    for j=1:J
        if abs(YY(j,1))<1e-6
            YY(j,1) = 0;
        end
    end
    %VV XX FF分别表示WIOT中的增加值、中间投入以及最终消费的格式。
    XXi = XX.*(repmat(YY,[1,J]))./(repmat(max(YY-NN,0),[1,J]));%repmat的意思是把YY按照1*J的维度平铺。
    XXi(isnan(XXi)|isinf(XXi)) = 0;
    
    FFi = zeros(J,N*4);
    for n=1:N
        FFi(:,(((n-1)*4))+1:(n*4)) = FF(:,(((n-1)*5)+1):(((n-1)*5)+4));
    end
    FFi = FFi.*(repmat(YY,[1,N*4]))./(repmat(max(YY-NN,0),[1,N*4]));
    FFi(isnan(FFi)|isinf(FFi)) = 0;
    FFdi = zeros(J,N); 
    for n=1:N
        FFdi(:,n) = sum(FFi(:,(((n-1)*4)+1):(n*4)),2);%去掉了存货
    end
      FFdi = FFdi;
     
      k = find(FFdi<0)
       FFdi(FFdi<0)=0;
    VVi = zeros(J,1);%增加值
    for n1 = 1:J
        VVi(n1) = sum(XXi(n1,:)) + sum(FFdi(n1,:)) - sum(XXi(:,n1));%说明sum中的为向量时，计算为所有向量的和。若为向量，且指定维度为1或2时，则按照列和按行加总。
    end
    VVi(VVi<0)=0;

    %IO Coefficients and Value Added Coeff | do net inventory corrections
    gammas = []; 
    for n= 1:N
        for s = 1:S
            inputlist = XXi(:,(n-1)*S+s);
            inputlist = reshape(inputlist,S,N);
            inputlist = sum(inputlist,2);
            inputlist = [inputlist;VVi((n-1)*S+s)];
            inputratio = inputlist./sum(inputlist);   
            gammas = [gammas,inputratio];
        end
    end
    VAcoeff = gammas(S+1,:); 
    VAcoeff(isnan(VAcoeff)) = 1;
    VAcoeff = reshape(VAcoeff,S,N);
    gammas = gammas(1:(size(gammas,1)-1),1:(size(gammas,2))); %size(a,b)b为1即为行数，b为2即为列数。
    % Problematic columns: Have all gammas=NaN (when output equals to 0);
    % Assume VAcoeff=1 for these industries, and gammas =0;
    gammas(isnan(gammas)) = 0;
    
    %XXi(XXi==0) = 1e-18;
    %FFdi(FFdi==0) = 1e-18;
     
    %Consumption Shares - Alphas
    alphas = zeros(S,N);
    CtyCons = sum(FFdi,1)'; %41*1，这个国家消费的所有最终品
    for n=1:N
        for s = 1:S
            alphas(s,n) = sum(FFdi(s:S:(N-1)*S+s,n))/CtyCons(n);
        end
    end

    % Shares
    Piagg = zeros(S,J);
    PiFagg = zeros(S,N);
    
    %20200515新加程序，为了计算Ij，第一种简单方法。
    XXi = XXi.*tauX;
    FFdi =FFdi.*tauF;%sum后为2意思是计算行向量的值。
    XXi_guanshui = sum(XXi.*(tauX-1))';%计算中间品关税
    FFdi_guanshui = sum((FFdi.*(tauF-1))')';%计算最终品关税
    XXi_gs= sum(reshape(XXi_guanshui,S,N))'; 
    FFdi_gs= sum(reshape(FFdi_guanshui,S,N))';    
    
   %下面不是新加
    for s=1:S
        Piagg(s,:) = sum(XXi([s:S:(N-1)*S+s],:),1);
        PiFagg(s,:) = sum(FFdi([s:S:(N-1)*S+s],:),1);
    end
    
%     for n1=1:N
%         Piagg = Piagg + XXi((S*(n1-1))+1:n1*S,:);
%         PiFagg = PiFagg + FFdi((S*(n1-1))+1:n1*S,:);
%     end
    Pi = XXi./repmat(Piagg,N,1);
    PiF = FFdi./repmat(PiFagg,N,1);
    Pi(isnan(Pi))=0;
    PiF(isnan(PiF))=0;
    F1_test=1-sum(PiF./tauF.*(tauF-1).*repmat(alphas,N,1));
    CtyConsp_F_test=CtyCons./F1_test';
    IGS21 = sum(Pi./tauX.*(tauX-1).*repmat(gammas,N,1));
    IGS22 = reshape(IGS21,56,44);
    IGS3=sum(IGS22.*reshape(YY,56,44));
    %Surplus and VA per country
    E = zeros(N,1);
    M = zeros(N,1);

    % for n1  = 1:1:N
    %     % Imports
    %     n1s =((n1-1)*S+1):n1*S;
    %     M(n1) = sum(sum(XXi(:,n1s))) - sum(sum(XXi(n1s,n1s))) + sum(FFdi(:,n1)) - sum(FFdi(n1s,n1));
    %     % Exports
    %     E(n1) = sum(sum(XXi(n1s,:))) - sum(sum(XXi(n1s,n1s))) + sum(sum(FFdi(n1s,:))) - sum(FFdi(n1s,n1));
    % end
    % I have to keep the Sn consistent - and this is the only way - since we
    % have used 
    %Imports
    M = zeros(N,1);
    %Exports
    E = zeros(N,1);
    E_test = zeros(N,1);
    M_test = zeros(N,1);
    for j = 1:N
        for s = 1:S
            for i = 1:N
                iiss = (i-1)*S+s;
                jjss = (j-1)*S+s;  
                for r = 1:S
                    jjrr = (j-1)*S+r;
                    iirr = (i-1)*S+r;
                    M(j) = M(j) + Pi(iiss,jjrr)/tauX(iiss,jjrr)*gammas(s,jjrr)*YY(jjrr);
                    E(j) = E(j) + Pi(jjss,iirr)/tauX(jjss,iirr)*gammas(s,iirr)*YY(iirr);
                    M_test(j) = M_test(j) + Pi(iiss,jjrr)/tauX(iiss,jjrr)*gammas(s,jjrr)*YY(jjrr);
                    E_test(j) = E_test(j) + Pi(jjss,iirr)/tauX(jjss,iirr)*gammas(s,iirr)*YY(iirr);
                end
                M(j) = M(j) + PiF(iiss,j)/tauF(iiss,j)*alphas(s,j)*CtyCons(j);
                E(j) = E(j) + PiF(jjss,i)/tauF(jjss,i)*alphas(s,i)*CtyCons(i);
                M_test(j) = M_test(j) + PiF(iiss,j)/tauF(iiss,j)*alphas(s,j)*CtyConsp_F_test(j)+PiF(iiss,j)*alphas(s,j)*IGS3(j)/F1_test(j);
                E_test(j) = E_test(j) + PiF(jjss,i)/tauF(jjss,i)*alphas(s,i)*CtyConsp_F_test(i)+PiF(jjss,i)*alphas(s,i)*IGS3(i)/F1_test(i);
            end
        end
    end
    %Surplusn = E-M;
    Surplusn=E_test-M_test;
  
    VAn = zeros(N,1);
    for n = 1:N
        VAn(n) = sum(VVi(((n-1)*S+1):n*S));
    end
    
    %Verify that Surplusn = VAn - CtyCons (Balanced)
end 
