function [ys, params, check]=tank_steadystate(ys, exo, M_, options_)

global  betac sigmac sigman chik epsilon_n deltac_d deltac_c alphac_c alphac_d phi_k phi_y epsilon_y epsilon_yy ...
        mu_dmg delta_x rho_r vi_pi vi_y rho_c rho_d ac_ss ad_ss r_ss pi_ss

number_of_parameters = M_.param_nbr;

for ii=1:number_of_parameters
    paramname=M_.param_names{ii};
    eval([ paramname ' = M_.params(' int2str(ii) ');']);
end
check = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%
options=optimset('MaxFunEvals',10000,'TolX',1e-10,'TolFun',1e-10);
Mh=betac;
%r=1/betac;
rkd=(1/betac)+deltac_d-1;
rkc=(1/betac)+deltac_c-1;
%xd=1;
%xd1=1;
%qd=1;
%xc=1;
%xc1=1;
%qc=1;
%mc=(epsilon_yy-1)/epsilon_yy;
%pi=1;
ac=1; 

kd_g =12.0201; %% 8.10; 
kc_g = 5.42985; %%5.5;;
nd_g=0.3; %0.277195; %0.8;
nc_g=0.3; %%0.173317;  %%0.4 %%0.6;
%%em_g=0.483104; %%0.5;
cl_g=0.7;
 
guess = [kd_g, kc_g, nd_g, nc_g, cl_g]; %%cl_g
ss_sol = fsolve(@findss_tank, guess, options);


kd = ss_sol(1);
kc = ss_sol(2);
nd=ss_sol(3);
nc=ss_sol(4);
cl=ss_sol(5);

%%n=(nc^(1+epsilon_n)+nd^(1+epsilon_n))^(1/(1+epsilon_n));
id=deltac_d*kd;
ic=deltac_c*kc;
ec=(kc^alphac_c)*(nc^(1-alphac_c));  %%%damage on this?
ed=(kd^alphac_d)*(nd^(1-alphac_d));

%x=ed/(1-delta_x);
%dmg=1-exp(-mu_dmg*x);
y=ac*(ec^((epsilon_y-1)/epsilon_y)+ed^((epsilon_y-1)/epsilon_y))^(epsilon_y/(epsilon_y-1));  %%(1-dmg)*
%pc = (y/ec)^(1/epsilon_y);
%pd = (y/ed)^(1/epsilon_y);
ch=y-cl-id-ic;
c=cl+ch;
wd=(1-alphac_d)*(ed/nd);
wc=(1-alphac_c)*(ec/nc);
gini_income=(wc*nc)/(wd*nd);
gini_consumption=(ch/cl);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

params=NaN(number_of_parameters,1);
for  iter = 1:length(M_.params) %update parameters set in the file
     eval([ 'params(' num2str(iter) ') = ' M_.param_names{iter} ';' ])
end

number_of_endogenous_variables = M_.orig_endo_nbr;
for ii = 1:number_of_endogenous_variables
    varname = M_.endo_names{ii};
    eval(['ys(' int2str(ii) ') = ' varname ';']);
end