function [T_wat,H_lake,q]=water_only_new(T_wat,SWR_in,cc,T_a,U,V,h_alt,RH,H_lake)

del_t=86400;
rho_w=1000;
% rho_i=917;
rho_i=rho_w;
J=1.907*10^(-5);
Lf=3.34*10^5;
Cw=4186;
Kw=0.56;
Ki=2.24;
sigma=5.67*10^(-8);

%   Newton Raphson module

Ts=T_wat(1);
N=T_wat(2);
for k=1:30

    T_v=[Ts-1 Ts Ts+1];
    for j=1:3
        
        I0=0.6;
        
        %   Energy balance of the lake
        
        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))));   %   Saturation vapour pressure of air
        e=0.01*RH*e_s;                                      %   Vapour pressure of air
        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)))));  %   Saturation vapour pressure at surface
        e_surf=e_s_s*0.01*RH;                                   %   Vapour pressure at surface        
        H=rho_a*(1000)*Ct*u*(T_a-T_v(j));
        LE=rho_a*(2.501*10^6)*Ct*u*(e-e_surf)*(0.622/P);        
        if u==0
            H=0;
            LE=0;
        end
        Fc=1000*4184*J*sign(N-T_v(j))*(abs(N-T_v(j)))^(4/3);
        Q(j)=SWR_in*(1-I0)+LWR_in-sigma*T_v(j)^4+H+LE+Fc;
        if j==2
            q=SWR_in*(1-I0)+LWR_in-sigma*(T_v(j))^4+H+LE;
        end
    end
    tsn=Ts-Q(2)/((Q(3)-Q(1))/2);
    Ts=tsn;
%    disp(Ts);
end
T_wat(1)=Ts;
end
