function [inhead,inlake,staghead1,Nt6,TSRt6,SSt6,SRt6,hhSt6,hhRt6,Vl_t,Qcl_t]=filling_flood_mod4U_tM(ohsrcfx,srfslpgrd,xaxs,yaxs,Qlk,wc100,t13,levar,Vt_obs);
%[inhead,inlake,staghead,Vl_t,Qcl_t]=filling_flood_mod4U_t(ohsrcfx,srfslpgrd,xaxs,yaxs,lklist,wc100,t13,levar,hots,plotzon);
%by Sasha Carter 10-17-13
%allows variable input into your lake based ont the output of the ssflux
%models
%updated Sept 1, 2015
%now uses filling_flood_mod6UH where I had a better treatment of Alpha,
%this also means that I need to update levar and make sure I have it right
%additionally I have improved the conservation of mass the sheetflow and
%possibly added sediment advection
km=1e3;
[~,nc]=size(Qlk);
Lnumcol=3;%Lake number
Ql4=4;%the extra stuff that's part of Qlake
LQcol=5;%influx into lake from upstream
ISAT_t=[3.835259882;4.193939386;4.406192041;4.790116527;5.190235634;5.525538311;5.988653284;6.297628722;...]
    6.534492873;6.846306116;7.230529054;7.808150303;8.184288546;8.771460363;8.930059528;9.266375648;9.749705134;
    10;10.5;11;11.5;12;12.5;13;13.5;14.1;14.7;15.3;15.8;16.2];
[ncIS]=length(ISAT_t);%VALIDATION
fI0=find(ISAT_t>=t13(1)&ISAT_t<=t13(2));%
ISAT_t0=ISAT_t(fI0);
[ncIS0]=length(ISAT_t0);%VALIDATION_0
ISAT_t1=[t13(1);ISAT_t0;t13(2)];
[ncIS1]=length(ISAT_t1);%VALIDATION_1
QI0a=max(fI0(1)-1,1);%which Q intervals are you going to use?
QI0=[QI0a;fI0];%
if t13(2)>max(ISAT_t)
    QI0(length(QI0))=QI0(length(QI0)-1);
%else
%    QI0=[QI0a;fI0];%
end
t13a=[ISAT_t1(1:2)];
Qlk4now=Qlk(:,[1:Ql4,Ql4+QI0(1)]);
hotfac=0;
pltfac=0;
c=1;
[inhead0,inlake0,staghead0,Nt6,TSRt6,SSt6,SRt6,hhSt6,hhRt6,Vl_t,Qcl_t,c]=filling_flood_mod6UH(ohsrcfx,srfslpgrd,xaxs,yaxs,Qlk4now,wc100,t13a,levar,[],pltfac,c);
disp(inlake0)
disp(Qlk4now)
flk4inl=find(inlake0(:,Lnumcol)==levar(8));
fQlk4n=find(Qlk4now(:,Lnumcol)==levar(8));
[nrV,ncV]=size(Vl_t);
Vl_t(:,ncV)=[];
Vrg=max(Vl_t(2:nrV,:),[],2)-min(Vl_t(2:nrV,:),[],2);
fV4n=find(Vrg==max(Vrg));
figure(5)
plot(Vt_obs(1,:),Vt_obs(2,:),'g',Vl_t(1,:)+2e3,Vl_t(1+fV4n,:)/(km^3),'b');
xlabel('year');
ylabel('volume km^3');
legend('OBS','MOD');
disp('(Vl_t1)')
disp(size(Vl_t))
overhead52
disp(size(Vl_t));
for n1=2:ncIS1-1
    figure(4)
    plot(-ohsrcfx(:,1)/km,ohsrcfx(:,6));
     t13a=[ISAT_t1(n1:n1+1)];
     %disp(QI0)
     disp([flk4inl,Ql4+1,size(inlake0),size(Qlk),fQlk4n,Ql4+QI0(n1)]);
     inlake0(flk4inl,Ql4+1)=Qlk(fQlk4n,Ql4+QI0(n1));
     hotfac=1;
     %pltfac=0;
     [inhead1,inlake1,staghead1,Ntt6,TSRtt6,SStt6,SRtt6,hhStt6,hhRtt6,Vl_t1,Qcl_t1,c]=filling_flood_mod6UH(inhead0,srfslpgrd,xaxs,yaxs,inlake0,wc100,t13a,levar,staghead0,pltfac,c);
    %[inhead,inlake,staghead,Nt,Vl_t,Qcl_t,wcrow]=filling_flood_mod6UH(ohsrcfx,srfslpgrd,xaxs,yaxs,lklist,wc100,t13,levar,staghead0,plotzon,wcrow);
     [~,ncV0]=size(Vl_t1);
     [~,ncQ0]=size(Qcl_t1);
     [~,ncvlt]=size(Vl_t);
     Vl_t=[Vl_t(:,1:ncvlt-1),Vl_t1(:,2:ncV0-1)];
     disp(size(Vl_t));
     Qcl_t=[Qcl_t,Qcl_t1(:,1:ncQ0-1)];
     inhead0=inhead1;
     staghead0=staghead1;
     inlake0=inlake1;
     figure(5)
     plot(Vt_obs(1,:),Vt_obs(2,:),'g',Vl_t(1,:)+2e3,Vl_t(1+fV4n,:)/(km^3),'b');
     xlabel('year');
     ylabel('volume km^3');
     legend('OBS','MOD');
     overhead52
     Nt6=[Nt6,Ntt6];%later this preserves the option of more detailed
     TSRt6=[TSRt6,TSRtt6];%later this preserves the option of more detailed
     SSt6=[SSt6,SStt6];%later this preserves the option of more detailed
     SRt6=[SRt6,SRtt6];%later this preserves the option of more detailed
     hhRt6=[hhRt6,hhRtt6];
     hhSt6=[hhSt6,hhStt6];
     %output
     disp('(Vl_t2)')
    disp(size(Vl_t))
     %ned=input('hit enter now')
 end
 inhead=inhead0;
 inlake=inlake0;