function [del_I,T_s_i,H_lake,virt_lid_flag]=refreeze(T_wat,SWR_in,cc,T_a,U,V,h_alt,RH,H_lake,del_I,virt_lid_flag,scal,T_s_i)
     
     rho_w=1000;
     rho_i=rho_w;
     del_t=86400;
     sigma=5.67*10^(-8);
     J=1.907*10^(-5);
     Lf=3.34*10^5;
     Cw=4186;
     Ki=2.24;
     melt=0;
     
     Ts=T_a;

     for i=1:25
	     T_v=[Ts-1 Ts Ts+1];
	     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*RH*0.01;
		     H=rho_a*(1000)*Ct*u*(T_a-T_v(j));
		     LE=rho_a*(2.835*10^(6))*Ct*u*(0.622/P)*(e-e_surf);
		     if u==0
			     H=0;
			     LE=0;
		     end
		     GHR=-(Ki/(H_lake/(scal)+del_I))*(-T_wat(2)+T_v(j));
		     Q(j)=SWR_in+LWR_in-sigma*((T_v(j))^4)+H+LE+GHR;
             Q1(j)=SWR_in+LWR_in-sigma*((T_v(j))^4)+H+LE;
	     end
	     tsn=Ts-Q(2)/((Q(3)-Q(1))/2);
	     Ts=tsn;
     end

     kdTdz=Ki*(-T_wat(2)+Ts)/((H_lake/(scal))+del_I);
     NN=-(kdTdz*del_t)/(rho_i*Lf);
     T_s_i=Ts;
%      disp(Ts);
    if Ts<273.16    % Further freezing
        if NN>0
            del_I=del_I+NN;
            if NN<H_lake
                H_lake=H_lake-NN;
            else
                H_lake=0;
            end
        end
        virt_lid_flag=1;
    else   %    Melting takes place
        T_s_i=273.16;
        kdTdz=Ki*(T_s_i-273.16)/((H_lake/(scal))+del_I);
        M=0;
        Q_f=freezing_temp(SWR_in,cc,T_a,U,V,h_alt,RH);
        NN=((Q_f-kdTdz)*(del_t))/(rho_i*Lf);
        if NN>0
            if NN>del_I  %  Whole lid melts
                H_lake=H_lake+del_I;
                del_I=0;
            else         %  Part of the lid melts
                H_lake=H_lake+NN;
                del_I=del_I-NN;
            end
        end
        if del_I<=0
            virt_lid_flag=0;
            del_I=0;
        end
    end
end
