function [T_s_i,T_wat,H_lake,H_ice]=ice_water_new_3(H_lake,U,V,SWR,cc,RH,T_a,T_s_i,T_wat,h_alt,H_ice,k1,scal)
del_t=86400;
alpha_s=0.431;
Ls=2.835*10^(6);
Lf=3.34*10^(5);
Ki=2.24;
rho_w=1000;
sigma=5.67*10^(-8);
Cpa=1005;
ci=2100;
Cw=4180;
J=1.907*10^(-5);
rho_i=rho_w;
Q1=0;
Melt=0;


%   Estimating melt at bottom of the ice-lid and at the lake base

    if T_wat(2)>273.16 && H_lake>=0.1
        T1=-((del_t*Ki)/(rho_i*Lf))*((273.16-T_wat(2))/(H_lake/scal));
        T2=sign(-273.16+T_wat(2))*((Cw*rho_w*del_t*J)/(rho_i*Lf))*(abs(-273.16+T_wat(2)))^(4/3);
%         T2=0;
        del_h_bel=T1+T2;
%         H_lake=H_lake+del_h_bel;
    else
        T1=((del_t*Ki)/(rho_i*Lf))*((T_s_i-273.16)/(H_ice+H_lake/scal));
        T2=0;
        del_h_bel=T1+T2;
%         del_h_bel=0;
    end
%     T_s_i=273.16;
    if T_wat(2)>273.16 && H_lake>=0.1
        
        T1=((del_t*Ki)/(rho_i*Lf))*((T_s_i-273.16)/(H_ice+H_lake/scal));
%         T1=((del_t*Ki)/(rho_i*Lf))*((T_wat(2)-273.16)/(H_lake/scal));
%         T1=((del_t*Ki)/(rho_i*Lf))*((T_s_i-T_wat(2))/(H_ice+H_lake/scal));
        T2=sign(-273.16+T_wat(2))*((Cw*rho_w*del_t*J)/(rho_i*Lf))*(abs(-273.16+T_wat(2)))^(4/3);
    else
       T1=((del_t*Ki)/(rho_i*Lf))*((T_s_i-273.16)/(H_ice+H_lake/scal));
       T2=0;
    end

    del_h_ice_bel=T1+T2;        %   Ice-lid base

    H_ice=H_ice-del_h_ice_bel;
    H_lake=H_lake+del_h_ice_bel+del_h_bel;
    
if T_wat(2)>273.16 && H_lake>=0.1

%   Estimation of water layer temperatures
    
    A=rho_w*Cw*sign(T_wat(2)-273.16)*(abs(T_wat(2)-273.16))^(4/3);
    B=rho_w*Cw*sign(T_wat(2)-273.16)*(abs(T_wat(2)-273.16))^(4/3);
    T_wat(2)=T_wat(2)-(del_t*J/(rho_w*Cw*H_lake))*(A+B);
else
    T_s_i=T_s_i-(Cw/ci)*(H_lake/H_ice)*(273.16-T_wat(2));
    T_wat(2)=273.16;
end

%   Newton Raphson loop

Ts=T_s_i;
for k=1:25
    
    %   Estimation of ice surface temperature
    
    T_v=[Ts-1 Ts Ts+1]; %   Assume no melt first
    for j=1:3
        u=sqrt(U^2+V^2);
        rb=(9.8*(T_a-T_v(j))*10)/(T_a*(u^2));
        if rb<0
            Ct=(0.0013)*(1-(2*20*rb)/(1+50.986*sqrt(abs(rb))));
        else
            Ct=(0.0013)*(1+20*rb)^(-2);
        end       
        P=101325*exp(-8000/h_alt);
        rho_a=(P*0.02896)/(8.31*T_a);
        e_s=611*exp((2494/0.461)*((1/273.16)-(1/(T_a))));
        e=0.01*RH*e_s;
        LWR_in=((0.23+0.484*(e/T_a)^(1/8))*(1-cc^4)+0.952*cc^4)*sigma*(T_a^4);
        e_s_s=611*exp((2494/0.461)*((1/273.16)-(1/(T_v(j)))));
        e_surf=e_s_s*0.01*RH;
        H=rho_a*Cpa*Ct*u*(T_a-T_v(j));
        LE=rho_a*Ls*Ct*u*(e-e_surf)*(0.622/P);
        if u==0
            H=0;
            LE=0;
        end

        GHR=((Ki)/H_ice/(1))*(273.16-T_v(j));
        Q(j)=SWR*(1-alpha_s)+LWR_in-sigma*(T_v(j))^4+H+LE+GHR;
    end
    tsn=Ts-Q(2)/((Q(3)-Q(1))/2);
    Ts=tsn;
   
end

if Ts>=273.16

    T_s_i=273.16;
    u=sqrt(U^2+V^2);
    rb=(9.8*(T_a-273.16)*10)/(T_a*(u^2));
    if rb<0
        Ct=(0.0013)*(1-(2*20*rb)/(1+50.986*sqrt(abs(rb))));
    else
        Ct=(0.0013)*(1+20*rb)^(-2);
    end
    P=101325*exp(-8000/h_alt);
    rho_a=(P*0.02896)/(8.31*T_a);
    e_s=611*exp((2494/0.461)*((1/273.16)-(1/(T_a))));
    e=0.01*RH*e_s;
    LWR_in=((0.23+0.484*(e/T_a)^(1/8))*(1-cc^4)+0.952*cc^4)*sigma*(T_a^4);
    e_s_s=611*exp((2494/0.461)*((1/273.16)-(1/273.16)));
    e_surf=e_s_s*0.01*RH;
    H=rho_a*Cpa*Ct*u*(T_a-273.16);
    LE=rho_a*Ls*Ct*u*(e-e_surf)*(0.622/P);
    if u==0
        H=0;
        LE=0;
    end
    GHR=((Ki)/H_ice)*(273.16-273.16);
    Q1=SWR*(1-alpha_s)+LWR_in-sigma*(273.16)^4+H+LE+GHR;
    Melt=(Q1*del_t)/(rho_i*Lf);
    H_lake=H_lake+Melt;
    H_ice=H_ice-Melt;
else
    T_s_i=Ts;
end

end
