function [inhead,inlake,staghead,Nt,TSRt,SSt,SRt,hhSt,hhRt,Vl_t,Qcl_t,wcrow]=filling_flood_mod6UH(ohsrcfx,srfslpgrd,xaxs,yaxs,lklist,wc100,t13,levar,staghead0,plotzon,wcrow);
% [inhead,inlake,Vl_t,Qcl_t]=filling_flood_mod1(ohsrcfx_SLM2_6,srfslp_1kf,xax1k,yax1k,Q8_09K1_2d,wc100,0,3.5);
%levar is [Qfloodon,Qfloodoff,v_s,D_s,N_inf,R_kl,ntstep]
%[inhSLWs2,ilkSLWs2,shSLW1,VlSLWs1,QSLWs1]=filling_flood_mod4(ohsrcfx_SLW_3,srfslp_1kf,xax1k,yax1k,Q8_06K1_2d,wc100,0,5);
%staghead,dheadds,dheaddt,ustag,ustag_old
%staghead,dheadds,dheaddt,ustag,Srf_t,tvout,smod
%this program give you a sheetflow model that eventually leads to an
%outburst
%ohsrcfx is the file of the form
%[s,x,y,bed,phead,head,lrep,w3,w5,w10,drno,row,col,melt,lkno,Qb,sideflux]
%This version is based on the work of Kinglsake and Ng 2013 Ann glaciol
%includes both a sheet and a Rothlisberger channel
%starts withthe output of your ssflux algorithim ideally should not be too
%different but sometimes ends up being so
%latest uupdate:  uses some walder fowler sediment erosion concepts to
%allow the conduit to grow via mechanical erosion and contract via sediment
%deformation, sediment deformeation assumed to be identical to ice except
%of a different vales for n and for K
%tries an upwinding algorithim to get flux and sheet dimensions and
%ultimately 
%plotzon=3; %shows every few steps% this is now in the inputs
ttplot=1; %days between plots
plot_fon=1;%plots final results 
%inhead=inhead0;
[nr0,nc0]=size(ohsrcfx);
%define obvious paramters
scol=1; %distance downstream
xcol=2;%x(m)
ycol=3;%(y(meters))
bedcol=4;%bed elevation (m)
phdcol=5;
hcol=6;%hydro head
lrcol=7;
w3col=8;
w5col=9;
brkcol=10;
strncol=11;
lkcol=12;
rowcol=13;
colcol=14;
bcol=15;%meltrate
Qbcol=16;%baseflow
dQbcol=17;%in and out Baseflow
%END OF WHAT COMES FROM OHSRFCX
tkcol=18;%ice thickness
srfcol=19;%ice surface
Scol=20;%aperature
hscol=21;%sheet thickness
Ncol=22;%effective pressure (Pas)

pocol=23;%overburden pressure (Pas)
pwcol=24;%water pressure (Pas)
hpascol=25;%hdopotential in pascals
Eshtcol=26;%
Tshtcol=27;

Vscol=28;%sheetflow velocity :)
frcol=29;

edgcol=30;%new as of 2013 = 1 if a lake edge 0 otherwisealways inside the lake boundary
fstatcol=31;%is it filling 0 is not a lake, 1 filling, 2 draining, 3 streamflood
storcol=32;%storage along this section of flow
Taubcol=33;%basal shear stress
gsslpcol=34;%grid cell slope
flslpcol=35; %surface slope along flow
cellidcol=36;
mc2col=37;%old melt rate
Tsrcol=38;%distributed - channel transmission
Rccol=39;%conduit radius
Tauccol=40;%Conduit stress for erosion
Dccol=41;%deposition
Krpccol=42;%creep closure rate
Eccol=43;%Erosion
Qccol=44;%conduit flux
Nccol=45;%conduit effective pressure
Sccol=46;%conduit X-sectional area
Cbarccol=47;%Conduit sed concentration;
Vccol=48;%conduit avg water velocity
pwccol=49;%conduit water pressure
hpasccol=50;%conduit hydropotential
Kccol=51;%hydraulic radius
htestccol=52;% conduit reynolds number
mccol=53;%conduit wall melt
Recol=54;

ncf=55;
conducol=Qccol:Recol;
stgcol=[scol:Tsrcol,Nccol,hpasccol,htestccol,Recol];
format short g
% this is where we start with the original output of the model
% upstream end is 1st line hopefully a lake is there
% hopefully a lake or the ocean is at the bottom
%now define the lake columns
Lxcol=1;
Lycol=2;
Lnumcol=3;%Lake number
Lareacol=4;%aerial extent
LQcol=5;%influx into lake from upstream
Lhavcol=6;%average hydrootential of lake 
LRcol=7;%lake radius--important for a really crude flexure I am trying to do 
LCcol=8;%s position of the lake center
Lsealcol=9;% position of the cell with the highest hydropotential downstream of the lake
Luscol=10;%upstream most cell
Ldscol=11;%downstreammost cell 
Lhpmincol=12;%lowest apparent hydropotential of lake
Lvolcol=13;%Lake volume
nclkf=13;

%discharge cms
%define constants
tsat0=t13(1);
tstop0=t13(2);
nyr=1/52;
dayyr=365; %days per year
secday=3600*24;%seconds per day
secyr=secday*dayyr;
Ai=6.8e-24;% ps3s(1/3)  Glenns flow law A
Bi=5.28e7;%pas(1/3) Glenns flow law B
Bclap=-7.4e-8;%K/pa Cassius claperon Beta coefficient
Bcomp=1e-7;%copressibility gof water Pa^-1
Nprime=0.1; % manning roughness coefficient of Tmp=to/Bclap*pw
Aroth=7.27e7; %rothlisbergers exponent
Glen_n=3; %glenns flow law exponents
Rau_w=1e3;%density of water kg/m3
Rau_i=917; % denisty of ice
Rau_s=2.7e3; %density of sediment
g=9.81; %gravitational acceleration
cp_w=4217.7; %J/(Kg*K)
Lh_fus=333500;%J/kg heat fo fusion
Kh_topo=0.0309; %topo correction
Kw=0.558;% Thermal conductivity of water W/mK

Ecap=0.313;%pattersons melt number
pmtc=-8.7e-4;%pressure melting point depression degrees K per meter
pmpp=9.8e-8;%same as above but in K/pa
mu_h2o=1.787e-3;%Pas-1
Prdt=mu_h2o*cp_w/Kw;%prandtl number
K1=(2/27)*(4*pi())^(2/3);
rothrat=11/(8*Glen_n);%Magic ratio used all over the place;
nprime=0.1;
Banana=(1/Bi)*(1/nprime)^(6/11)*(Glen_n*Aroth)^(8*Glen_n/11);

%from kingslake et al., 2013;
Rkkl=0.005;%transfer efficiency between conduits and channels
kkl=1e-9;
K1kl=1.1;
pkl=4;
qkl=1;
cskl=2e-20;%m s-1 Pa-3
%fkl=0.07%orig
fkl=0.07;%m-2/3s2  THIs one may vary around too
K0kl=1e-24*1;%pa-3s-1
K0klc=1e-24*1;%pa-3s-1
K0klf=1e-24*300;%pa-3s-1
Rkl=levar(6);
%Rkl=0.012;%this is going to really vary around a bit  FUDGEABLE
%was 0.006
%now onto the walder and fowler constants 1994
WFa=1.33;
WFb=1.8;
WFAs=3e-5;%sediment deformability
WFfr=0.1;%friction factor
WFK1=.1;
WFK2=6;
%WFvs=levar(3);%settling velocity based on an assumed median grain sie
WFvs=levar(4)^2*(2/9)*(Rau_s-Rau_w)*g/(mu_h2o*levar(11));%this is the better way to do it
%WFalpha=levar(4)^2*(2/9)*(Rau_s-Rau_w)*g/(levar(3)*mu_h2o)%algebraicly
%reworked
%WFvs=5e-6;
%WFDs=1.2e-4;%characteristic grain size
WFDs=levar(4);
%WFD50=2.5e-4;%median grain size
WFD50=2*WFDs;%median grain size
%neededfor Kl&ng  stolen from Clarke 2003 and Flowers 200X
RR_not=.2;%initial radius
mfdg=100;
%0.012 is from PIG trunk (Rippin et al., 2011 J glaciol)
%0.1 is from Walkder 1986
km=1e3;
N_inf=levar(5)*Rau_w*g;%*effective pressure of till
%N_inf=g*Rau_w*9;%*effective pressure of till  right now at 10 m w.e  
Tau_ca=0.00125*WFD50*g*(Rau_s-Rau_i)*1; %critical shear stress
%Tau_c=0.208;
Tau_c=0.0052;
disp('Tau_c')
disp([Tau_c,Tau_ca]);
Q_floodon=levar(1);%m3s-1Coefficient at which your flood becomes real 
Q_floodoff=levar(2);%m3s-1Coefficient at which your flood becomes real
Q_proto=-0.5;%initial flood outflow
LGEOfac=levar(9);

Rhonnie=nthroot((2/pi)^5,3)*Nprime*sqrt(Rau_w*g)/2;
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];
[ncIS]=length(ISAT_t);%VALIDATION


%now some operational stuff  YOU CAN ADJUST THESE 04/29/13
tsfil=secyr/52; %month
tsdrain=secday;
ntstep=levar(7);
ntrec=25;
ntdisp=150;%how frequently do you want the cfl timestep and other stuff displaying on the screen
fntsntr=floor(ntstep/ntrec);
ntCFL=20;
%Banana=(1/Bi)*(nprime)^(6/11)*(Glen_n*Aroth)^(8*Glen_n/11);
tstep=secday/24;%timestep = 1 hour
nstep=floor(nyr*dayyr*secday/tstep);%how long the model runs
k_n=5e-4;
%k_n=5e-3;
grsz=2.5e3;%2-d grid size
ohssz=1e3;%1-d step size
%SO FAR IT DOES NOT MATTER
%k_n=0;
Nmax=10*Rau_w*g;%convert from m.w.e. to pascals.
lkofkey=levar(8);
disp(levar)
%Now we start to maniupulate this mess
%********BEGIN HERE***********
if isempty(staghead0)==1
    %disp('coldstart')
    %ned=input('hit enter for me');
    [nr0,nc0]=size(ohsrcfx);
    smod=[ohsrcfx(:,1),zeros(nr0,ncIS)];%more VALIDATION
    [nrl0,ncl0]=size(lklist);
    
    nrl2=0;
    lklist3a=rmrept(sortrows(ohsrcfx(:,lkcol)));%News as of 05/2013 lakelist improved
    lklist3=lklist3a(lklist3a>0.5);%sorts all the different lakes in lakecol
    %until you have a list complere with no repeat
    disp(lklist3);
    nlake=length(lklist3);
    lklist2=zeros(nlake,nclkf);
    %05/13 makeing inlake Re did this so that every lake makes it on to lakelist
    %the ice shelf cavity at present should not rise as water flow into it
    %we can acomplish this by making it ridiculously big or settign NET inflow to zero
    for n1=1:nlake;
        lk4now=lklist3(n1);
        disp(lk4now)
        flk=find(abs(lklist(:,Lnumcol)-lk4now)<=0.1);
        %if isempty(flkohsrc)==0&&lk4now>1
        flkohsrc=find(abs(ohsrcfx(:,lkcol)-lk4now)<=0.1);
        mnlkdex=min(flkohsrc);
        mxlkdex=max(flkohsrc);
        lakebit=zeros(1,nclkf);
        lakebit(Luscol)=mxlkdex;%height of highets point in lake
        lakebit(Ldscol)=mnlkdex;%average height of lake
        if isempty(flk)==0
            lakebit(Lxcol:LQcol)=lklist(flk,Lxcol:LQcol);
            if mnlkdex==1;
                lk_diam=ohsrcfx(mxlkdex,scol)-ohsrcfx(mnlkdex,scol);
                lakebit(Lsealcol)=1;
                lakebit(LRcol)=lk_diam/2;
                lakebit(LCcol)=(ohsrcfx(mnlkdex,scol)+ohsrcfx(mxlkdex,scol))/2;
                lakebit(Lsealcol)=1;
                lakebit(Lhpmincol)=min(ohsrcfx(mnlkdex:mxlkdex,hcol))*Rau_w*g;
            else
                
                lk_diam=ohsrcfx(mxlkdex,scol)-ohsrcfx(mnlkdex-1,scol);
                lakebit(LRcol)=lk_diam/2;
                lakebit(LCcol)=median(ohsrcfx(mnlkdex:mxlkdex,scol));
                disp(mnlkdex)
                lkbelowdex=(find(ohsrcfx(1:mnlkdex-1,lkcol)>0.5,1,'last'));
                sealh4now=max(ohsrcfx(lkbelowdex+1:mnlkdex-1,hcol));
                %disp(ohsrcfx(1:mnlkdex-1,lkcol)>0.5)
                lakebit(Lsealcol)=find(ohsrcfx(1:mnlkdex-1,hcol)==sealh4now);
                lakebit(Lhpmincol)=min(ohsrcfx(mnlkdex:mxlkdex,hcol))*Rau_w*g;
                %disp(lakebit)
                
            end
            
        else%somehow the lake is not on our list WhOOPS!
            %lakebit(Lxcol:LQcol)=lklist(flk,Lxcol:LQcol);
            lakebit(Lxcol:Lnumcol)=[median(ohsrcfx(flkohsrc,xcol:ycol),1),lk4now];
            lakebit(Lareacol)=mean(ohsrcfx(flkohsrc,w3col))*abs(ohsrcfx(mnlkdex,scol)-...
                ohsrcfx(mxlkdex,scol));%aerial extent
            lakebit(LQcol)=abs(ohsrcfx(mxlkdex,Qbcol));%make the in
            if mnlkdex==1;
                lk_diam=ohsrcfx(mxlkdex,scol)-ohsrcfx(mnlkdex,scol);
                lakebit(Lsealcol)=1;
                lakebit(LRcol)=lk_diam/2;
                lakebit(LCcol)=(ohsrcfx(mnlkdex,scol)+ohsrcfx(mxlkdex,scol))/2;
                lakebit(Lsealcol)=1;
                lakebit(Lhpmincol)=min(ohsrcfx(mnlkdex:mxlkdex,hcol))*Rau_w*g;
                
            else
                lk_diam=ohsrcfx(mxlkdex,scol)-ohsrcfx(mnlkdex-1,scol);
                lakebit(LRcol)=lk_diam/2;
                lakebit(LCcol)=median(ohsrcfx(mnlkdex:mxlkdex,scol));
                disp(mnlkdex)
                lkbelowdex=(find(ohsrcfx(1:mnlkdex-1,lkcol)>0.5,1,'last'));
                sealh4now=max(ohsrcfx(lkbelowdex+1:mnlkdex-1,hcol));
                %disp(ohsrcfx(1:mnlkdex-1,lkcol)>0.5)
                lakebit(Lsealcol)=find(ohsrcfx(1:mnlkdex-1,hcol)==sealh4now);
                lakebit(Lhpmincol)=min(ohsrcfx(mnlkdex:mxlkdex,hcol))*Rau_w*g;
            end
            nrl2=nrl2+1;
            
        end
        lakebit(Lhavcol)=lakebit(Lhpmincol);
        if abs(lk4now-1)<.1||floor(lk4now/1111)==lk4now/1111%are we under an ice shelf?
            lakebit(LQcol)=0;%the ocean is removing water as fast as it enters so net influx = 0
        end
        lklist2(n1,:)=lakebit;
    end
    disp(lklist2)
    
    
    %begin working with inhead
    inhead=[ohsrcfx,zeros(nr0,ncf-nc0)];
    inlake=[lklist2];
    keylake0=inlake(inlake(:,Lnumcol)==lkofkey,:);
    lkdsdex0=find(inhead(1:keylake0(Lsealcol),hcol)<keylake0(Lhpmincol)/(Rau_w*g),1,'last');
    lkdsdex2=find(abs(ohsrcfx(:,lkcol)-keylake0(Lnumcol))<0.25,1,'first');
    dscell=[lkdsdex0:lkdsdex2];%do I use this?
    dscellsg=[lkdsdex0:lkdsdex2-1];%staggered lake is smaller for rchannel
    fdestlake0=find(ohsrcfx(1:lkdsdex2-1,lkcol)>0.5,1,'last');
    destlakenum=ohsrcfx(fdestlake0,lkcol);
    rachel=[fdestlake0+1:max(dscellsg)];%the portion fot he flowpath below the main lake but above the next lake
    rachel0=min(rachel);
    rachelf=max(rachel);
    rachels0=rachel0-1;
    rachelsg=[fdestlake0:max(dscellsg)];%the staggered grid extends the reach into the lake
    %identifying the lake with largest inflow and then theupstream most point
    %with lower head so that I can know when there is a monotonically
    %decreasign path between the lake and that point
    %xy2_s(ohsrcfx(:,2:nc0),1,2);
    %make sure we have decent
    %taub b is importatn here and we decided to filter it over  teh scale of a
    %ffew ice thicknesses as otherwise it gets kind of messy
    headmat=zeros(nr0*2-1,ncf);
    gslp0=gttrack(srfslpgrd,xaxs,yaxs,ohsrcfx(:,xcol:ycol));
    gsfltr0=triang(7);
    gsfltr=gsfltr0/sum(gsfltr0);
    gslp=filtfilt(gsfltr,[1],gslp0);%we need a smoother surface slope
    inhead(:,gsslpcol)=gslp(:,3);
    
    for i0=1:nr0;
        headmat(2*i0-1,:)=inhead(i0,:);
    end
    [~,pkrow]=findpeaks(ohsrcfx(:,Qbcol),'MINPEAKDISTANCE',5);% this distance = the resolution of the melt model = 5 km (joughin et al., 2004)
    sqbpre=ohsrcfx([1;pkrow;nr0],[scol,Qbcol]);
    sqbpost=pchip(sqbpre(:,1),sqbpre(:,2),ohsrcfx(:,scol));
    
    inhead(:,scol)=-ohsrcfx(:,scol);%necessary becasue kinglake and ng 2013 do istance backwards realtive to
    inhead(:,tkcol)=inhead(:,phdcol)*Rau_w/Rau_i; %ice thickness
    inhead(:,srfcol)=inhead(:,tkcol)+inhead(:,bedcol);%surface elevation
    inhead(:,pocol)=inhead(:,phdcol)*Rau_w*g;% overburden pressure
    inhead(:,hpascol)=inhead(:,hcol)*Rau_w*g;%hydropotential
    inhead(:,Taubcol)=(inhead(:,pocol).*inhead(:,gsslpcol));%basal shear stress
    inhead(:,cellidcol)=inhead(:,rowcol)*1e5+inhead(:,colcol);
    inhead(2:nr0-1,edgcol)=real((inhead(2:nr0-1,lkcol)>=.5).*(inhead(1:nr0-2,lkcol)<=.5|inhead(3:nr0,lkcol)<=.5));%lake edges parameterized:)
    %disp(inhead(:,[lkcol,edgcol]))
    if inhead(2,lkcol)<0.5&&inhead(1,lkcol)>0.5
        inhead(1,edgcol)=1;
    end
    if inhead(nr0-1,lkcol)<0.5&&inhead(nr0,lkcol)>0.5
        inhead(nr0,edgcol)=1;
    end
    
    %new material to account for side inflow and melt input kinda
    inhead(:,Qbcol)=sqbpost;%m3/s;
    %inhead(:,dQbcol)=[sqbpost(2:nr0)-sqbpost(1:nr0-1);inhead(nr0,Qbcol)]/ohssz;%m2/s oe m3/sper km length  divide by km to get per m unit length
    justds=[(inhead(2:nr0,scol)-inhead(1:nr0-1,scol))];%on the staggered grid
    disp(size(justds));
    CFLdx=min(justds);
    justdsg=[(inhead(2,scol)-inhead(1,scol));(inhead(3:nr0,scol)-inhead(1:nr0-2,scol))/2;(inhead(nr0,scol)-inhead(nr0-1,scol))];%on the regular grid
    
    vwide=grsz*ones(nr0,1);%when in doubt 2.5 km is the width
    %while satis<2 %debugging now commented out
    dhdskl0=[inhead(2,:)-inhead(1,:);((inhead(3:nr0,:)-inhead(1:nr0-2,:)));inhead(nr0,:)-inhead(nr0-1,:)];
    
    for i2=1:nr0;
        ds4now31=inhead(min([nr0,i2+1]),scol)-inhead(max([i2-1,1]),scol);
        inhead(i2,flslpcol)=((inhead(min([nr0,i2+1]),srfcol)-inhead(max([i2-1,1]),srfcol))/ds4now31);
        %dhds=((inhead(i2+1,hcol)-inhead(i2-1,hcol))/ds4now31);
        dhdskl=dhdskl0(i2,:)/ds4now31;
        vwide(i2)=min([grsz*4,max([inhead(i2,w3col),2e2])]);%tdid it this way because of 5 km grid size
        %disp([vwide,dhds])
        %inhead(i2,hscol)=abs(nthroot(inhead(i2,Qbcol)*12*mu_h2o/(Rau_w*g*dhds*vwide(i2)),3));
        
        %disp([i2,inhead(i2,[hscol,Ncol])])
        %inhead(i2,Ncol)=min([Nmax,k_n*inhead(i2,Taubcol)/(Rau_w*g*abs(max([inhead(i2,hscol),5e-5])))]);
        lk4now=inhead(i2,lkcol);
        
        if lk4now<=0.5%storage within the cell ignores tortuosity
            %inhead(i2,Ncol)=min([Nmax,k_n*inhead(i2,Taubcol)/(abs(max([inhead(i2,hscol),5e-5])))]);
            inhead(i2,Scol)=abs((pi()*Rkl/(4*K1kl))^(-2/3)*sqrt(Rau_w*g)*sqrt(fkl/6.6)*inhead(i2,Qbcol)...
                *(dhdskl(hpascol))^(-1/2));%equation 9 from Kingslake and Ng 2013
            %threw an ab sign on there as negative channel geometries are nto
            %realisitc int his universe
            %/vwide(i2);
            %inhead(i2,hscol)=max(inhead(i2,Scol)/vwide(i2),ohsrcfx(i2,bedcol)-ohsrcfx(i2,brkcol));
            %this is the source of the instability at SLM.  It is worse
            %with modded hydropotential
            inhead(i2,hscol)=inhead(i2,Scol)/vwide(i2);
            inhead(i2,storcol)=inhead(i2,hscol)*vwide(i2)*abs(ds4now31/2);%
            inhead(i2,Ncol)=(((pi()*Rkl*cskl*Glen_n^Glen_n)/(4*K0kl))*(inhead(i2,Taubcol).^pkl)/inhead(i2,Scol))^(1/(Glen_n+qkl));
            %equation 11a in kingslake and Ng 2013
            inhead(i2,pwcol)=inhead(i2,pocol)-inhead(i2,Ncol);
            inhead(i2,hpascol)=(inhead(i2,bedcol))*Rau_w*g+inhead(i2,pwcol);
            inhead(i2,Scol)=inhead(i2,hscol)*vwide(i2);
            inhead(i2,fstatcol)=0;
            inhead(i2,brkcol)=inhead(i2,bedcol)-inhead(i2,hscol);
            %keep the ice base (bed) consist wiht measurments brk is by inference shouldn't matter much, but I want it to be this way
        else%treat "lakes differentlyu from non lake cells"  allow cells to change
            filake=find(abs(inlake(:,Lnumcol)-lk4now)<0.01);
            if inhead(i2,edgcol)==1
                inhead(i2,Scol)=abs((pi()*Rkl/(4*K1kl))^(-2/3)*sqrt(Rau_w*g)*sqrt(fkl/6.6)*inhead(i2,Qbcol)...
                    *(dhdskl(hpascol))^(-1/2));%equation 9 from Kingslake and Ng 2013
                inhead(i2,hscol)=inhead(i2,Scol)/vwide(i2);
                %inhead(i2,Ncol)=(((pi()*Rkl*cskl*Glen_n^Glen_n)/(4*K0kl))*(inhead(i2,Taubcol).^pkl)/inhead(i2,Scol))^(1/(Glen_n+qkl));
                inhead(i2,brkcol)=inhead(i2,bedcol)-inhead(i2,hscol);
                %equation 11a in kingslake and Ng 2013
            else
                inhead(i2,brkcol)=min([ohsrcfx(i2,brkcol),inhead(i2,bedcol)-0.1]);%ensures a nonzero lakedepth everywhere to eliminate instabilites
                %realted to zero water layer thicknesses
                inhead(i2,hscol)=inhead(i2,bedcol)-inhead(i2,brkcol);
            end
            if isempty(filake)==0%is it on the lake list?
                %lake is at least 1 m deep
                %disp(i2)
                inhead(i2,storcol)=lklist(filake,4)*1*abs(ds4now31/2);
                inhead(i2,Ncol)=ohsrcfx(i2,hcol)*Rau_w*g-inlake(filake,Lhpmincol);%adapted from 3b to make
                %lake flexmore  This apparently didn't work so hot
                %disp(inlake(filake,Lhpmincol)/(Rau_w*g))
                inhead(i2,pwcol)=inhead(i2,pocol)-inhead(i2,Ncol);
                inhead(i2,hpascol)=(inhead(i2,brkcol)+inhead(i2,hscol))*Rau_w*g+inhead(i2,pwcol);
                
                inhead(i2,fstatcol)=1;
            else%generally reserved for ocean cavities but maybe secondary lakes as well why not
                % disp(inhead(i2,[lkcol,edgcol]))
                %disp(i2)
                lkrow=find(inhead(:,lkcol)==lk4now);
                flkmin=find(inhead(lkrow,hpascol)==min(inhead(lkrow,hpascol)));
                %if i2~=lkrow(flkmin)
                %    inhead(i2,lkcol)=0;
                %end
                inhead(i2,storcol)=inhead(i2,hscol)*vwide(i2)*abs(ds4now31/2);%
                hpmin=min(inhead(lkrow,hcol));
                inhead(i2,Ncol)=ohsrcfx(i2,hcol)*Rau_w*g-hpmin*Rau_w*g;%adapted from 3b to make
                inhead(i2,pwcol)=inhead(i2,pocol)-inhead(i2,Ncol);
                inhead(i2,hpascol)=(inhead(i2,brkcol)+inhead(i2,hscol))*Rau_w*g+inhead(i2,pwcol);
                %disp(inhead(i2,hpascol))
                inhead(i2,fstatcol)=0;%through flowing lakes need some parameterization
            end
            if inhead(i2,edgcol)==0
                inhead(i2,Scol)=inhead(i2,hscol)*vwide(i2);
            end
            
            
        end
        %for now
        
    end
    lkrow0=find(ohsrcfx(:,lkcol)>=1);
    nolkrow=find(abs(ohsrcfx(:,lkcol))<0.01);
%     nolkrowc=nolkrow(nolkrow>1.5&nolkrow<(nr0-.5));
%     stnolkrow0=sort([nolkrow(nolkrow<(nr0-.5));inlake(:,Luscol)]);
%     vwide_sh=(vwide(1:nr0-1)+vwide(2:nr0))/2;
%     
    %need to make sure the fluxes like melt are distributed properly  %this is
    %a substantial change from previous model 4-30-13%
    %we smoothed the plot fo baselfo vs duatance and just take the chaneg in
    %smoothed baseflow along the path as the input/output
    inhead(:,dQbcol)=-[inhead(2:nr0,Qbcol)-inhead(1:nr0-1,Qbcol);0]/ohssz;%m2/s oe m3/sper km length  divide by km to get per m unit length
    %THIS IS WHERE YOU LEFT OFF MARCH 16, 2011
    inhead_preloop=inhead;
    %Begin Cutpaste here
    %all the materail below is there to record change and plot it
    inhead=inhead_preloop;
    staghead=zeros(nr0-1,ncf);%staggered grid on which we try to calcualte fluxes
    dheadds=zeros(nr0-1,ncf);%
    dheaddt=zeros(nr0,ncf);
    
    Vl_t=zeros(1+nlake,fntsntr+1);%volume of sourcelake
    %disp('vl_t')
    %disp(size(Vl_t))
    Qcl_t=zeros(3,fntsntr+1);%outflow from sourcelake
    %TSRt6,SSt6,SRt6,hhSt6,hhRt6
    Nt=zeros(1+nr0,fntsntr+1);
    TSRt=zeros(1+nr0,fntsntr+1);
    SSt=zeros(1+nr0,fntsntr+1);
    SRt=zeros(1+nr0,fntsntr+1);
    hhSt=zeros(1+nr0,fntsntr+1);
    hhRt=zeros(1+nr0,fntsntr+1);
    Qcl_t(1,1)=tsat0;
    Vl_t(1,1)=tsat0;
    %TSRt6,SSt6,SRt6,hhSt6,hhRt6
    Nt(1,1)=[tsat0];
    TSRt(1,1)=[tsat0];
    SSt(1,1)=[tsat0];
    SRt(1,1)=[tsat0];
    hhRt(1,1)=[tsat0];
    hhSt(1,1)=[tsat0];
    stag_lake=(inhead(1:nr0-1,lkcol)+inhead(2:nr0,lkcol))/2;
    CFLudex0=find(stag_lake==0);%the non lake cells
    CFLudex1=find(rem(stag_lake*2,2)==1);%the lake edges;
    %
    CFLudex2=sort([CFLudex0;CFLudex1]);
    CFLudex=CFLudex2;
    CFLudex(CFLudex2==nr0-1)=[];%kill the ends
    CFLudex(CFLudex2==1)=[];%kill the ends
    figure(7);
    plot(inhead(:,scol)/km,inhead(:,srfcol),inhead(lkrow0,scol)/km,inhead(lkrow0,srfcol),'b*');
   % figure(4);
    %plot(inhead(:,scol)/km,inhead(:,hpascol)/(Rau_w*g),inhead(lkrow0,scol)/km,inhead(lkrow0,hpascol)/(Rau_w*g),'b*');
    for i3=1:nr0-1
       staghead(i3,stgcol)=(inhead(i3,stgcol)+inhead(i3+1,stgcol))/2;
    end
    %
    nlake=length(lklist3);
    conduon=0;
    %plot_lake_geom(ohsrcfx);
    troncrit=0;
else
    %disp('hotstart')
    %ned=input('hit enter for me');
    inhead=ohsrcfx;
    inlake=lklist;
    staghead=staghead0;
    keylake0=inlake(abs(inlake(:,Lnumcol)-lkofkey)<.5,:);
    %
    justds=[(inhead(2:nr0,scol)-inhead(1:nr0-1,scol))];%on the staggered grid
    disp(size(justds));
    CFLdx=min(justds);
    justdsg=[(inhead(2,scol)-inhead(1,scol));(inhead(3:nr0,scol)-inhead(1:nr0-2,scol))/2;(inhead(nr0,scol)-inhead(nr0-1,scol))];
    %
    %keylake0=inlake(inlake(:,Lnumcol)==max(inlake(:,LQcol)),:);
    lkdsdex0=find(inhead(1:keylake0(Lsealcol),hcol)<keylake0(Lhpmincol)/(Rau_w*g),1,'last');
    lkdsdex2=find(abs(ohsrcfx(:,lkcol)-keylake0(Lnumcol))<0.25,1,'first');
    dscell=[lkdsdex0:lkdsdex2];%do I use this?
    dscellsg=[lkdsdex0:lkdsdex2-1];%staggered lake is smaller for rchannel
    fdestlake0=find(ohsrcfx(1:lkdsdex2-1,lkcol)>0.5,1,'last');
    destlakenum=ohsrcfx(fdestlake0,lkcol);
    rachel=[fdestlake0+1:max(dscellsg)];%the portion fot he flowpath below the main lake but above the next lake
    rachel0=min(rachel);
    rachelf=max(rachel);
    rachels0=rachel0-1;
    rachelsg=[fdestlake0:max(dscellsg)];%the staggered grid extends the reach into the lake
    %
    lkrow0=find(ohsrcfx(:,lkcol)>=1);
    nolkrow=find(abs(ohsrcfx(:,lkcol))<0.01);
    %
    vwide=grsz*ones(nr0,1);
    for i2=1:nr0;
        %dhds=((inhead(i2+1,hcol)-inhead(i2-1,hcol))/ds4now31);
        vwide(i2)=min([grsz*4,max([inhead(i2,w3col),2e2])]);
    end
    %
    lklist3a=rmrept(sortrows(ohsrcfx(:,lkcol)));%News as of 05/2013 lakelist improved
    lklist3=lklist3a(lklist3a>0.5)
    nlake=length(lklist3);
    %
    if max(abs(staghead(:,Qccol)))>Q_floodoff
       if abs(staghead(rachelf,Qccol))<Q_floodoff
           wtf83=max(abs(staghead(:,Qccol)));
           disp('wtf83')
           disp([wtf83,find(abs(staghead(:,Qccol))==wtf83)]);
       end
        conduon=1;
        troncrit=-2;
    else
        conduon=0;
        troncrit=0;
    end
    disp('conduon=')
    disp(conduon)
    %
    %staghead=zeros(nr0-1,ncf);%staggered grid on which we try to calcualte fluxes
    dheadds=zeros(nr0-1,ncf);%
    dheaddt=zeros(nr0,ncf);
    %
    Vl_t=zeros(1+nlake,fntsntr+1);%volume of sourcelake
    %disp('vl_t2')
    %disp(size(Vl_t))
    Qcl_t=zeros(3,fntsntr+1);%outflow from sourcelake
    %TSRt6,SSt6,SRt6,hhSt6,hhRt6
    Nt=zeros(1+nr0,fntsntr+1);
    TSRt=zeros(1+nr0,fntsntr+1);
    SSt=zeros(1+nr0,fntsntr+1);
    SRt=zeros(1+nr0,fntsntr+1);
    hhSt=zeros(1+nr0,fntsntr+1);
    hhRt=zeros(1+nr0,fntsntr+1);
    %Qcl_t(1,1)=tsat0;
    %Vl_t(1,1)=tsat0;
    Qstag_hot=staghead(:,Scol)*(pi()*Rkl/(4*K1kl))^(2/3)*sqrt(6.6/(Rau_w*g*fkl)).*dheadds(:,hpascol);
    Qcl_t(:,1)=[tsat0;Qstag_hot(rachelf-1);staghead(rachelf,Qccol)];
    Vl_t(:,1)=[tsat0;inlake(:,Lvolcol)];
    %disp('vl_t3')
    %disp(size(Vl_t))
    %TSRt6,SSt6,SRt6,hhSt6,hhRt6
    Nt(:,1)=[tsat0;inhead(:,Ncol)];
    TSRt(:,1)=[tsat0;inhead(:,Tsrcol)];
    SSt(:,1)=[tsat0;inhead(:,Scol)];
    SRt(:,1)=[tsat0;inhead(:,Sccol)];
    hhSt(:,1)=[tsat0;inhead(:,hpascol)];
    hhRt(:,1)=[tsat0;inhead(:,hpasccol)];
    
end
mo_justds=mode(justds)
%Bedrock 1
%inhead(:,brkcol)=inhead(:,bedcol);
trec1=1;
dstagdt=zeros(size(staghead));
%THIS IS WHERE THE MAIN LOOP BEGINS!
stopdex=0;
t4now=0;%seconds
dlakedt=zeros(size(inlake)); 
%wcrow=c;%plotting infor
wcrmax=length(wc100(:,1));

amifirst=1;
pigs=1;
fly=2;
stpplot=0;
dispon=0;

holyshit=0;%a trigger to halt the implosion of the model
%disp(keylake0(Lnumcol))

%ned=input('this is keylake')
for i4=2:ncf
    dheadds(:,i4)=(inhead(2:nr0,i4)-inhead(1:nr0-1,i4))./(inhead(2:nr0,scol)-inhead(1:nr0-1,scol));
end

%troncrit=2;
for j1=1:ntstep
    for i3=1:nr0-1
        
        staghead(i3,stgcol)=(inhead(i3,stgcol)+inhead(i3+1,stgcol))/2;
    end
    for i4=2:ncf
        dheadds(:,i4)=(inhead(2:nr0,i4)-inhead(1:nr0-1,i4))./justds;
    end
    if abs(sum(sign(dheadds(dscellsg,hpascol))))<length(dheadds(dscellsg,hpascol))&&abs(troncrit)<=1
        troncrit=0;%reset
    end
    if abs(sum(sign(dheadds(dscellsg,hpascol))))==length(dheadds(dscellsg,hpascol))&&abs(troncrit)<=1;%is there a monotonically decreasing path in the hpot to the next natural low
      %  disp('commence leakage of main lake')
      %  disp([j1,staghead(rachelf-1,Eshtcol)])
        troncrit=troncrit+0.05;%if you have enough in a row you get a flood
        
    end
    %if troncrit>.75
    %   break
    % end
    %this addresses what happens at the lake edge as the lake
    %2 ways to do this, either extend the channel into the lake and give it
    %continuous properties like a hose
    flkedge=find(floor(staghead(:,lkcol))~=(staghead(:,lkcol)));
    flkedge1=find(floor(staghead(:,lkcol))~=(inhead(1:nr0-1,lkcol)));
    flkedge2=find(floor(staghead(:,lkcol))~=(inhead(2:nr0,lkcol)));
    flkedge=rmrept(sort([flkedge1;flkedge2]));
    L_flkedge=length(flkedge);
    Qstag=staghead(:,Scol)*(pi()*Rkl/(4*K1kl))^(2/3)*sqrt(6.6/(Rau_w*g*fkl)).*dheadds(:,hpascol);%sheetflux on the staggered grid
    staghead(:,Vscol)=Qstag./staghead(:,Scol);%priming up for WF_model assumes straight path
    staghead(rachelsg,Tshtcol)=abs((1/8)*Rau_w*WFfr*staghead(rachelsg,Vscol));
    staghead(rachelsg,Eshtcol)=WFK1*WFvs*(max((staghead(rachelsg,Tshtcol)-Tau_c),0)/(g*(Rau_s-Rau_w)*WFDs)).^(3/2);
    %begin conduits HERE  Calcualte
    %if troncrit>=1&&max(staghead(rachelsg,Eshtcol))>1e-9;%spin up the conduits since they can't actually function if all empty
    
    if troncrit>=1&&abs(Qstag(rachelf-1))>Q_floodon;%spin up the conduits since they can't actually function if all empty
        disp('there she blows')
        dheadds_2lak=(inhead(rachels0,hpascol)-inhead(rachelf,hpascol))/(inhead(rachels0,scol)-inhead(rachelf,scol));
        disp('dheadds')
        disp(dheadds_2lak)
        disp('hs')
        disp(inhead(rachelf-1,hscol))
        
        disp('years passed')
        disp(t4now/secyr);
        disp([j1,staghead(rachelf-1,Eshtcol)])
    
        inhead(lkrow0,Qccol)=0;%no conduits in the lake; possible however at lake outlet
        inhead(lkrow0,Sccol)=0;
        inhead(lkrow0,mccol)=0;%no conduits therefore no conduiit melt; possible however at lake outlet
        inhead(:,Nccol)=inhead(:,Ncol);%for continuity and this will keep water from transferring
        inhead(:,hpasccol)=inhead(:,hpascol);
        dheadds(:,hpasccol)=dheadds(:,hpascol);
        staghead(:,Nccol)=staghead(:,Ncol);%for continuity and this will keep water from transferring
        staghead(:,hpasccol)=staghead(:,hpascol);%again some cumbersome steps at spinup to get initial conditioons right
        %now begin working on "the reach"
        dstghpas=-([0;(staghead(2:nr0-1,hpasccol)-staghead(1:nr0-2,hpasccol))/km;0]);
        lakebit=keylake0;
        lxit=max([lakebit(Ldscol)-1,lakebit(Lsealcol)]);
        %identify the staggered cells
        disp('Qstag(lxit)');
        disp(Qstag(lxit));
        %inhead(rachel,Qccol)=-0.5;
        %staghead(rachelsg,Qccol)=-0.5;
        %inhead(rachel,Qccol)=Qstag(lxit);
        inhead(rachel,Qccol)=Q_proto;
        %staghead(rachelsg,Qccol)=Qstag(lxit);
        staghead(rachelsg,Qccol)=Q_proto;
        inhead(rachel,Sccol)=1*(abs((fkl*Rau_w*g*inhead(rachel,Qccol).^2)./(dstghpas(rachel)))).^(3/8);%redo of eq 3 in kl&ng
        staghead(rachelsg,Sccol)=1*(abs((fkl*Rau_w*g*staghead(rachelsg,Qccol).^2)./(dheadds(rachelsg,hpascol)))).^(3/8);
      
        %so we want to start with a guess for discharge at the lake exit
        %based on the formulation here and a constant channele gometry in
        %later timesteps the channel geometry downstream will evolve, but
        %we will still be making this initial guess at each timestep if we
        %follow kingslake and ng
        %staghead(rachels0:rachelf,Qccol)=sign(dheadds(rachels0:rachelf,hpasccol)).*sqrt(abs((dheadds(rachels0:rachelf,hpasccol).*staghead(rachels0:rachelf,Sccol).^(8/3))/(fkl*Rau_w*g)))*1;
        %inhead(rachel0:rachelf,Qccol)=sign(dstghpas(rachel0:rachelf)).*sqrt(abs((dstghpas(rachel0:rachelf).*inhead(rachel0:rachelf,Sccol).^(8/3))/(fkl*Rau_w*g)))*1;
        staghead(rachels0:rachelf-1,Qccol)=sign(dheadds(rachels0:rachelf-1,hpasccol)).*sqrt(abs((dheadds(rachels0:rachelf-1,hpasccol).*staghead(rachels0:rachelf-1,Sccol).^(8/3))/(fkl*Rau_w*g)))*1;
        inhead(rachel0:rachelf-1,Qccol)=sign(dstghpas(rachel0:rachelf-1)).*sqrt(abs((dstghpas(rachel0:rachelf-1).*inhead(rachel0:rachelf-1,Sccol).^(8/3))/(fkl*Rau_w*g)))*1;
        staghead(rachelsg,Vccol)=staghead(rachelsg,Qccol)./staghead(rachelsg,Sccol);%NEW WF
        inhead(rachel,Vccol)=inhead(rachel,Qccol)./inhead(rachel,Sccol);%NEW WF
        
        %troncrit=-2;%now we turn off the primer switch
        conduon=1;%and turn on the conduit math
        amifirst=1;
        disp('j1     dtCFLmax      minraduis        t4now        Qnetout   Qstag(just downstream)  Qccol  melt_c)');
        %break
    end
    
    if conduon>=1%Much of this needs to move to after the "STEP"
        staghead(rachelsg,mc2col)=mfdg*(((1-Kh_topo)*dheadds(rachelsg,hpasccol)+(Kh_topo)*Rau_w*g*dheadds(rachelsg,bedcol)...
            -(Kh_topo)*dheadds(rachelsg,Nccol)).*staghead(rachelsg,Qccol))/Lh_fus;% a bit concerende becasue
        %staghead(rachelsg,mccol)=100*(dheadds(rachelsg,hpasccol).*staghead(rachelsg,Qccol))/Lh_fus;%the
        %begin WF parameters for sed
        staghead(rachelsg,Vccol)=staghead(rachelsg,Qccol)./staghead(rachelsg,Sccol);
        staghead(rachelsg,Rccol)=sqrt(staghead(rachelsg,Sccol)/(pi()));
        staghead(rachelsg,Tauccol)=abs((1/8)*WFfr*Rau_w*staghead(rachelsg,Vccol));%stress on conduit walls
        staghead(rachelsg,Eccol)=WFK1*WFvs*(max((staghead(rachelsg,Tauccol)-Tau_c),0)/(g*(Rau_s-Rau_w)*WFDs)).^(3/2);
        staghead(rachelsg,Dccol)=WFK2*WFvs*staghead(rachelsg,Cbarccol).*sqrt(g*(Rau_s-Rau_w)*WFDs./staghead(rachelsg,Tauccol));
        staghead(rachelsg,mccol)=Rau_s*pi()*staghead(rachelsg,Rccol).*(staghead(rachelsg,Eccol)-staghead(rachelsg,Dccol));%now adapted for WF94 mode
        %for j5=rachelf:-1:rachels0
            %staghead(rachelsg,Cbarccol)=max(staghead(rachelsg+1,Cbarccol)+staghead(rachelsg,mccol)./(Rau_w*abs(staghead(rachelsg,Qccol))),0);%sediment advection term is relly really simplified here 
            staghead(rachelsg,Cbarccol)=max(staghead(rachelsg+1,Cbarccol)+staghead(rachelsg,mccol)./(Rau_w*abs(staghead(rachelsg,Qccol))),0);
            %assumes that the change in sediment concentration between
            %timesteps is not significantly large  and tries to account for
            %advection to some degree
            %also assumes that advection of seds is instantaneous
        %end
        %old version
        %this melt rate does not consider changes in pressure melting point May have to fix this ARRGH!
        inhead(rachel,Qccol)=[(staghead(rachels0:rachelf-1,Qccol)+staghead(rachel0:rachelf,Qccol))/2];%this neutralizes melt rates at minima and maxima
        inhead(rachel,Vccol)=inhead(rachel,Qccol)./inhead(rachel,Sccol);%newWF
        
        inhead(rachel,Rccol)=sqrt(inhead(rachel,Sccol)/(pi()));
        inhead(rachel,Tauccol)=abs((1/8)*WFfr*Rau_w*inhead(rachel,Vccol).^2);%stress on conduit walls
       % disp(inhead(69:79,[Tauccol,Vccol,Qccol]));
        %plot(inhead(:,1)/km,[inhead(:,hpasccol)/km,inhead(:,hpascol)/km]);
       % disp('[Sccol,Qccol,Vccol,Tauccol,Vscol]');
       % disp([inhead(rachelf,[Sccol,Qccol,Vccol,Tauccol]),staghead(rachelf-1,Vscol)])
       % ned=input('hit a button to end it all now');
        inhead(rachel,Eccol)=WFK1*WFvs*(max((inhead(rachel,Tauccol)-Tau_c),0)/(g*(Rau_s-Rau_w)*WFDs)).^(3/2);
        inhead(rachel,Dccol)=WFK2*WFvs*inhead(rachel,Cbarccol).*sqrt(g*(Rau_s-Rau_w)*WFDs./inhead(rachel,Tauccol));
        inhead(rachel,mccol)=Rau_s*pi()*inhead(rachel,Rccol).*(inhead(rachel,Eccol)-inhead(rachel,Dccol));%now adapted for WF94 mode
        %for j5=rachelf:-1:rachel0
            inhead(rachel,Cbarccol)=max(inhead(rachel+1,Cbarccol)+mo_justds*inhead(rachel,mccol)./(Rau_w*abs(inhead(rachel,Qccol))),0);%sediment advection term is relly really simplified here 
            %inhead(rachel,Cbarccol)=max(inhead(rachel+1,Cbarccol).*staghead(rachelsg+1,Qccol)+inhead(rachel,mccol)./(Rau_w*abs(inhead(rachel,Qccol))),0);
            %inhead(j5,Cbarccol)=max(inhead(j5+1,Cbarccol)+inhead(j5,mccol)/(Rau_w*abs(inhead(j5,Qccol))),0);
            %assumes that the change in sediment concentration between
            %timesteps is not significantly large  and tries to account for
            %advection to some degree
            %also assumes that advection of seds is instantaneous
        %end
        dhpasCstg=(staghead(rachel0:rachelf,hpasccol)-staghead(rachels0:rachelf-1,hpasccol))./(staghead(rachel0:rachelf,scol)-staghead(rachels0:rachelf-1,scol));
        dbedCstg=(staghead(rachel0:rachelf,bedcol)-staghead(rachels0:rachelf-1,bedcol))./(staghead(rachel0:rachelf,scol)-staghead(rachels0:rachelf-1,scol));
        dNCstg=(staghead(rachel0:rachelf,Nccol)-staghead(rachels0:rachelf-1,Nccol))./(staghead(rachel0:rachelf,scol)-staghead(rachels0:rachelf-1,scol));
        inhead(rachel,mc2col)=mfdg*((1-Kh_topo)*(dhpasCstg)+Kh_topo*Rau_w*g*dbedCstg+Kh_topo*dNCstg).*inhead(rachel,Qccol)/Lh_fus;
        %inhead(rachel,mccol)=100*((dhpasCstg).*inhead(rachel,Qccol))/Lh_fus;
        %inhead(rachel,mccol)=((dhpasCstg).*inhead(rachel,Qccol))/Lh_fus;
        %calcualte on staggered grid transfer to regular grid
        staghead(rachelsg,Krpccol)=sign(staghead(rachelsg,Nccol)).*(WFAs*staghead(rachelsg,Sccol).*...
        ((abs(staghead(rachelsg,Nccol))/WFa).^(WFa))/(2*N_inf^(WFb)));%eqn 3.6 from WF94
        inhead(rachel,Krpccol)=sign(inhead(rachel,Nccol)).*(WFAs*inhead(rachel,Sccol).*...
        ((abs(inhead(rachel,Nccol))/WFa).^(WFa))/(2*N_inf^(WFb)));%eqn 3.6 from WF94
        dstagdt(rachelsg,Sccol)=(staghead(rachelsg,mccol))/Rau_s-staghead(rachelsg,Krpccol);
        dheaddt(rachel,Sccol)=inhead(rachel,mccol)/Rau_s-inhead(rachel,Krpccol);
        Krpolds=K0klf*staghead(rachelsg,Sccol).*staghead(rachelsg,Nccol).^Glen_n;
        Krpold=K0klf*inhead(rachel,Sccol).*inhead(rachel,Nccol).^Glen_n;
        dheaddt(rachel,Sccol)=inhead(rachel,mccol)/Rau_s-inhead(rachel,Krpccol);%eq 1 from kingslake and ng 2013
%          dstagdt(rachelsg,Rccol);
%          dheaddt(rachel,Rccol);
%         dstagdt(rachelsg,Rccol);%WF94 model
%         dheaddt(rachel,Rccol);%WF94 model
        if isfucked(inhead)>0.5
            disp('dr')
            break
        end
        if isfucked(staghead)>0.5
            disp('stagdr')
            break
        end
        %figure(2)
        %plot(inhead(:,scol)/km,inhead(:,Sccol));
        %if troncrit>=1&&max(staghead(rachelsg,Eshtcol))>1e-9;
        if troncrit>=1&&abs(Qstag(rachelf-1))>Q_floodon
            troncrit=-2;%now we turn off the primer switch
            ntdisp=1000;%until we get these bugs worked out 
        end

    end
    %SHEETFLOW / Canals in Kl&NG 2013
    inhead(2:nr0-1,Qbcol)=(Qstag(2:nr0-1)+Qstag(1:nr0-2))/2;
    dQstag=[[0;(Qstag(2:nr0-1)-Qstag(1:nr0-2));0]./justdsg];
    dheaddt(nolkrow,Scol)=inhead(nolkrow,dQbcol)+dQstag(nolkrow)-inhead(nolkrow,Tsrcol);%Eqn 12 Kingslake and Ng 2013
    %dheaddt(nr0,Scol)=0;%this worked now furhter up in the code
    dheaddt(nolkrow,hscol)=dheaddt(nolkrow,Scol)./vwide(nolkrow);
    dheaddt(nolkrow,bedcol)=dheaddt(nolkrow,hscol);
    dheaddt(nolkrow,srfcol)=dheaddt(nolkrow,hscol);
    hned=find(inhead(:,hscol)<0);
    %Beging Timestepping HERE
    %LAKEDYNAMICS
    for i5=1:nlake%treating the lakes differently, yes
        lakebit=inlake(i5,:);
        lxit=max([lakebit(Ldscol)-1,lakebit(Lsealcol)]);
        l4nc=[lakebit(Ldscol):lakebit(Luscol)]';
        l4ncin=[(lakebit(Ldscol)+1):(lakebit(Luscol)-1)]';
        stgl4nc=[lakebit(Ldscol):(lakebit(Luscol)-1)]';%identify the staggered cells
        %if abs(lakebit(Lnumcol)-keylake0(Lnumcol))<0.1
        if abs(lakebit(Lnumcol)-lkofkey)<0.1
            Qnet=lakebit(LQcol)+Qstag(lxit)+staghead(rachelf,Qccol);
            Qnetout=Qnet;
        else
            if lakebit(Ldscol)>=2%rememebr flux downstream is negative in this reference frame
                Qnet=-((staghead(min([lakebit(Luscol)+1,length(staghead(:,1))]),Qccol)+Qstag(lakebit(Luscol)+1))-Qstag(lakebit(Ldscol)-2))*0;
                Qdest=0;
%                 if lakebit(Lnumcol)==71207111;
%                     Qdest=Qnet;
%                     Qdest2=Qstag(lakebit(Ldscol)-1);
%                 end
                
            else
                %you're in the ocean which should be at constgant head, at least presently,  if your eally want to knock yourself out you could include a tined modela here
                %that would be pretty cool actually, but not at the time of
                %this comment
                Qnet=0;
            end
        end
        %Qnet=lakebit(LQcol);
        inhead(l4nc,Qbcol)=-Qnet;%for some reason this has to be so part of the weirdness with tLh_fuse reversed s scale
        %Qstag(stgl4nc)=-inlake(i5,LQcol);
        Qstag(stgl4nc)=-Qnet;
        dlakedt(i5,Lvolcol)=Qnet;
        %dlakedt(i5,Lhpmincol)=2*Qnet/inlake(i5,Lareacol);%parabolic uplift
        %dlakedt(i5,Lhavcol)=2*Qnet/inlake(i5,Lareacol);
        %ul_scale=-((l4nc-mean(l4nc))*2/length(l4nc)).^2+1;%scaling parabolic uplift
        dlakedt(i5,Lhpmincol)=LGEOfac*Qnet/inlake(i5,Lareacol);%flat uplift
        dlakedt(i5,Lhavcol)=LGEOfac*Qnet/inlake(i5,Lareacol); 
        ul_scale=1;%scaling parabolic uplift
        %ulin_scale=ul_scale(2:(length(l4nc-1)));
        dheaddt(l4nc,hscol)=dlakedt(i5,Lhavcol)*ul_scale;%eq6 kl&ng'13
        dheaddt(l4nc,bedcol)=dheaddt(l4nc,hscol);%=hs
        dheaddt(l4nc,srfcol)=dheaddt(l4nc,hscol);%=hs
        dheaddt(l4nc,hpascol)=dlakedt(i5,Lhavcol)*(Rau_w*g);
        dheaddt(l4nc,hpasccol)=dlakedt(i5,Lhavcol)*(Rau_w*g);%conduit=sheet in lake
        dheaddt(l4nc,Ncol)=-dlakedt(i5,Lhavcol)*(Rau_w*g)*ul_scale;
        dheaddt(l4nc,Nccol)=-dlakedt(i5,Lhavcol)*(Rau_w*g)*ul_scale;%conduit=sheet in lake
        inhead(max(l4nc),Qbcol)=-Qnet;
    end
    fAsplode=find((abs(inhead(nolkrow,Sccol)./dheaddt(nolkrow,Sccol)))==nanmin(abs(inhead(nolkrow,Sccol)./dheaddt(nolkrow,Sccol))));
    if conduon>=1
        dtCFLmax_c=(1/20)*nanmin(abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol)));
        dtCFLmax_s=(1/20)*nanmin(abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol)));
        dtCFLmax=min([dtCFLmax_c,dtCFLmax_s,3600*6]);
        
    else
        dtCFLmax=(1/20)*min([nanmin(abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol))),3600*1]);
    end
    if dtCFLmax<2
        disp('implosion')
        if conduon<0.5
            dtCFLmax_s=dtCFLmax;
            dtCFLmax_c=dtCFLmax_s+1;
        end
        if dtCFLmax_c<dtCFLmax_s%if the conduit stops working then we transition back to sheets only
            fAsplode=find((abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol)))==nanmin(abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol))));
            disp(['conduit'])
            disp([fAsplode,rachel(fAsplode),inhead(rachel(fAsplode),[scol,Sccol]),dheaddt(rachel(fAsplode),Sccol)]);
            conduon=0;
            troncrit=0;
            staghead(:,[Qccol,Sccol,Nccol,Tsrcol,mccol])=0;
            inhead(:,[Qccol,Sccol,Nccol,Tsrcol,mccol])=0;
            dheaddt(:,Sccol)=0;
            dstagdt(rachelsg,Sccol)=0;
            
        else
            fAsplode=find((abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol)))==nanmin(abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol))));
            disp(['csheet']);
            disp([j1,nolkrow(fAsplode),inhead(nolkrow(fAsplode),[scol,Scol]),dheaddt(nolkrow(fAsplode),Scol)]);
        end
        %dispon=2;
        %stpplot=2;
    end
    if conduon>=1&&abs(inhead(rachelf-1,Qccol))<Q_floodoff%since the conduits turn on at a much higher Q this should be safe for now
            %fAsplode=find((abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol)))==nanmin(abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol))));
            disp(['drainage low  flood ends'])
            %disp([fAsplode,rachel(fAsplode),inhead(rachel(fAsplode),[scol,Sccol]),dheaddt(rachel(fAsplode),Sccol)]);
            conduon=0;
            troncrit=0;
            staghead(:,[Qccol,Sccol,Nccol,Tsrcol,mccol])=0;
            inhead(:,[Qccol,Sccol,Nccol,Tsrcol,mccol])=0;
            dheaddt(:,Sccol)=0;
            dstagdt(rachelsg,Sccol)=0;
    end
    if floor(j1/ntdisp)==(j1/ntdisp)%this is where we make plots of evolution over time so we can verify the model is behavign as planned
        %figure(7);
        %disp([(Qstag(rachelf-1)),Q_floodon])
        % line(inhead(:,scol)/km,inhead(:,srfcol),'color',wc100(wcrow,:));
        %figure(4);
        %line(inhead(:,scol)/km,inhead(:,hpascol)/(Rau_w*g),'color',wc100(wcrow,:));
        figure(5)
        line(t4now/secyr+t13(1)+2*km,inlake(abs(inlake(:,Lnumcol)-lkofkey)<.5,Lvolcol)/(1e9),'Marker','.');
        if wcrow>=wcrmax
            wcrow=1;
        else
            wcrow=wcrow+1;
        end
        
        if conduon>=1
% 
%                         figure(10)
%                         plot(inhead(:,scol)/km,[inhead(:,[hpasccol,hpascol])/(Rau_w*g),inhead(:,hcol)],...
%                             staghead(rachelsg,scol)/km,staghead(rachelsg,hpasccol)/(Rau_w*g)+staghead(rachelsg,Qccol),...
%                             staghead(rachelsg,scol)/km,-staghead(rachelsg,Qccol)+min(inhead(:,hcol)),...
%                             inhead(rachel,scol)/km,inhead(rachel,Tsrcol)*km+min(inhead(:,hcol)),...
%                             staghead(:,scol)/km,-Qstag+min(inhead(:,hcol)),...
%                             inhead(:,scol)/km,[min(inhead(:,hcol))+(inhead(:,mccol)*(1/Rau_w-1/Rau_i)+...
%                             K0kl*inhead(:,Sccol).*inhead(:,Nccol).^Glen_n)*1e5,min(inhead(:,hcol))*ones(nr0,1)]);
%                         title('hpot')
%                         legend('h_R','h_s','h_0','HpasC-Q_c','Q_c','Q_sht','Tsr','other dQ')
%                         axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,min(inhead(:,hcol))-5,max(inhead(rachel,hcol))+5]);
           % disp([j1,dtCFLmax,min(inhead(rachel,Sccol)),t4now/secyr,Qnetout,Qstag(rachelf-1),staghead(rachelf,Qccol),staghead(rachelf,mccol)*km,inhead(132:133,hscol)']);
           %disp('lo')
           disp([j1,dtCFLmax,t4now/secyr+t13(1),Qnetout,Qstag(rachelf-1),staghead(rachelf,Qccol),staghead(rachelf,mccol)*km,inhead(rachelf-1:rachelf,Sccol)',inhead(rachelf-1:rachelf,Nccol)']);
              %ned=input('push the button');
        else
            %disp([j1,dtCFLmax,min(inhead(:,Scol)),find(inhead(:,Scol)==min(inhead(:,Scol)),1,'first'),Qnetout,Qstag(rachelf-1),inhead(132:133,hscol)']);
            %fAsplode=find(abs((abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol)))-nanmin(abs(inhead(rachel,Sccol)./dheaddt(rachel,Sccol))))<1e-4,1,'first')
            %disp('yo')
            fAsplode=find((abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol)))==nanmin(abs(inhead(nolkrow,Scol)./dheaddt(nolkrow,Scol))));
            disp([fAsplode,j1,dtCFLmax,min(inhead(:,Scol)),t4now/secyr+t13(1),Qnetout,Qstag(rachelf-1),inhead(rachelf-1:rachelf,hscol)',inhead(rachelf-1:rachelf,Ncol)']);
        end
        %disp(hned);
    end
    if floor(j1/ntrec)==(j1/ntrec)%recording change over time
        trec1=trec1+1;
        Qcl_t(:,trec1)=[tsat0+t4now/secyr;Qstag(rachelf-1);staghead(rachelf,Qccol)];
        Vl_t(:,trec1)=[tsat0+t4now/secyr;inlake(:,Lvolcol)];
        %disp('vl_t4')
    %disp(size(Vl_t))
    %TSRt6,SSt6,SRt6,hhSt6,hhRt6
        Nt(:,trec1)=[tsat0+t4now/secyr;inhead(:,Ncol)];
        
    TSRt(:,trec1)=[tsat0+t4now/secyr;inhead(:,Tsrcol)];
    SSt(:,trec1)=[tsat0+t4now/secyr;inhead(:,Scol)];
    SRt(:,trec1)=[tsat0+t4now/secyr;inhead(:,Sccol)];
    hhSt(:,trec1)=[tsat0+t4now/secyr;inhead(:,hpascol)];
    hhRt(:,trec1)=[tsat0+t4now/secyr;inhead(:,hpasccol)];
    end
    inhead_old=inhead;
    if j1==1
        inhead_1=inhead;
    end
    %inhead(nolkrow,[Scol,hscol])=inhead(nolkrow,[Scol,hscol])+dheaddt(nolkrow,[Scol,hscol])*dtCFLmax;
    inhead=inhead+dheaddt*dtCFLmax;
    t4now=t4now+dtCFLmax;
    staghead(rachelsg,Sccol)=staghead(rachelsg,Sccol)+dstagdt(rachelsg,Sccol)*dtCFLmax;
    inlake=inlake+dlakedt*dtCFLmax;
    for i3=1:nr0-1
        
        staghead(i3,stgcol)=(inhead(i3,stgcol)+inhead(i3+1,stgcol))/2;
    end
    for i4=2:ncf%for some reason this needs to get calculated twice
        dheadds(:,i4)=(inhead(2:nr0,i4)-inhead(1:nr0-1,i4))./justds;
    end
    inhead(nolkrow,Ncol)=(((pi()*Rkl*cskl*Glen_n^Glen_n)/(4*K0kl))*(inhead(nolkrow,Taubcol).^pkl)./inhead(nolkrow,Scol)).^(1/(Glen_n+qkl));%Eq 11a Kl&ng13
    inhead(nolkrow,pwcol)=inhead(nolkrow,pocol)-inhead(nolkrow,Ncol);
    inhead(nolkrow,hpascol)=(inhead(nolkrow,brkcol)+inhead(nolkrow,hscol))*Rau_w*g+inhead(nolkrow,pwcol);
    
    if conduon>=1
        Qcprev2=staghead(rachelsg,Qccol);
        Nprev2=inhead(rachel,Nccol);
        dQcdsrec=zeros(length(rachelsg),1);
        dNcdsrec=dQcdsrec;
        %so hwere we go:  we calcualte N first, we have a guess for Q
        Nrmin=0;%effective pressure of a conduit cannot be less than zero
        hpotdest=inhead(rachels0,hpasccol);
        ffstep=1;
        %pkcrit=0;
        if abs(amifirst-1)<0.1;
            staghead(rachelf,Qccol)=staghead(rachelf,Qccol)*1;
            amifirst=0;
        else
            n4=rachelf;
            %dhpasc4now=(inhead(n4,hpasccol)-inhead(n4+1,hpasccol))/justds(n4);
            %Q_alt=-sign(dhpasc4now)*sqrt(abs((dhpasc4now*staghead(n4,Sccol)^(8/3))/(fkl*Rau_w*g)));
            Q_alt=Qguess;
            staghead(n4,Qccol)=Q_alt;
            Qguess=Q_alt;
        end
        Qguess=staghead(rachelf,Qccol);
        satis=0;
        itercount=0;
        ih_preiter=inhead;
        sh_preiter=staghead;
        %while satis<3
        for nsat=1:12% this is the slowest step and I wwish I could makeit more dynamic and adaptive
            Qcprev=staghead(rachelsg,Qccol);
            Nprev=inhead(rachel,Nccol);
            guesserr=zeros(1,3);
            guessstep=5e-4;
            for nn=1:3%begin the newton iteration
               % disp(nn)
                staghead(rachelf,Qccol)=Qguess+(nn-2)*guessstep;%guess+-dguess
                
                dhpasC_ideal=zeros(length(rachelsg),1);
                hpottydeal=zeros(length(rachelsg)*2,3);
                for n4=rachelf:-1:rachels0;
                    n41=n4-rachels0+1;
                    hpascmin=hpotdest;%the hydropotential should never be below that of the destination lake
                    dhpascmin=(hpascmin-inhead(n4+1,hpasccol))/(inhead(n4,scol)-inhead(n4+1,scol));
                    dhpascds0=fkl*Rau_w*g*(staghead(n4,Qccol)*abs(staghead(n4,Qccol)))/(staghead(n4,Sccol)^(8/3));%equation 3
                    dhpascds1=max([dhpascmin,dhpascds0]);%make sure that we don't go below destination lake
                    hpascmax=inhead(n4,hcol)*Rau_w*g-5e3;
                    dhpascmax=(hpascmax-inhead(n4+1,hpasccol))/(inhead(n4,scol)-inhead(n4+1,scol));
                    dhpascds=min([dhpascmax,dhpascds1]);
                    dhpascds=dhpascds0;
                    dNcdsrec(n4-rachels0+1)=dhpascds;%not sure why I have this might be vestigal
                    inhead(n4,hpasccol)=inhead(n4+1,hpasccol)+dhpascds*(inhead(n4,scol)-inhead(n4+1,scol));%calc hpot of next cell downstream
                    inhead(n4,Nccol)=(inhead(n4,hscol)+inhead(n4,brkcol))*Rau_w*g+inhead(n4,pocol)-inhead(n4,hpasccol);%calcualte effective pressure
                    inhead(n4,Tsrcol)=Rkkl*kkl*(inhead(n4,Nccol)-inhead(n4,Ncol));%calcualte transfger between conduit and sheet system
                    dhpasc4now=(inhead(n4,hpasccol)-inhead(n4+1,hpasccol))/(inhead(n4,scol)-inhead(n4+1,scol));
                    Qr_drive=staghead(n4,Qccol);
                    dQcds4=(inhead(n4,mccol)*(1/Rau_w-1/Rau_s)*0+inhead(n4,Krpccol)+inhead(n4,Tsrcol)/1);%eq 2  --Note the zero as mobillizing sediment doesn't change water volume
                    dQcdsrec(n4-rachels0+1)=dQcds4;
                    staghead(n4-1,Qccol)=staghead(n4,Qccol)+dQcds4*justdsg(n4-1);
                    dQ4now=(staghead(n4-1,Qccol)-staghead(n4,Qccol))/justdsg(n4-1);
                    N_old=inhead(n4,Nccol);
                    inhead(n4,hpasccol)=(inhead(n4,brkcol)+inhead(n4,hscol))*Rau_w*g+inhead(n4,pocol)-inhead(n4,Nccol);
                    inhead(n4,Tsrcol)=Rkkl*kkl*(inhead(n4,Nccol)-inhead(n4,Ncol));
                    %disp
                    dhpasC_ideal(n41)=fkl*Rau_w*g*staghead(n4,Qccol)*abs(staghead(n4,Qccol))/(staghead(n4,Sccol))^(8/3);%here we go with more fun and games
                    dhpasC_wegot=(inhead(n4+1,hpasccol)-inhead(n4,hpasccol))/(inhead(n4+1,scol)-inhead(n4,scol));
                    hpottydeal([(n41)*2-1,(n41)*2],1:3)=flipud([[inhead(n4+1,scol),[dhpasC_ideal(n41),dhpasC_wegot]*(inhead(n4+1,scol)-staghead(n4,scol))];...
                        [inhead(n4,scol),[dhpasC_ideal(n41),dhpasC_wegot]*(inhead(n4,scol)-staghead(n4,scol))]]);
                    if isfucked(inhead)>0.5
                        disp('bvp')
                        disp(itercount)
                        holyshit=1;
                        break
                    end
                    if isfucked(staghead)>0.5
                        disp('stagbvp')
                        disp(itercount)
                        holyshit=1;
                        break
                    end
                    %  ned=input('push a button to end it all now');
                end
                if holyshit==1
                    break
                end
                    if isfucked(inhead)>0.5
                        disp('bvp2')
                        disp(itercount)
                        holyshit=1;
                        break
                    end
                    if isfucked(staghead)>0.5
                        disp('stagbvp2')
                        disp(itercount)
                        holyshit=1;
                        break
                    end
                inhead(1:rachels0,htestccol)=inhead(rachels0,hpascol);
                for n5=rachels0:rachelf%integrate from destination lake to calcualte source lake hydropotential
                    n5a=n5-rachels0+1;
                    %disp([n5,n5a,2*n5a,2*n5a-1,hpottydeal(2*n5a,2),hpottydeal(2*n5a-1,2)])
                    inhead(n5+1,htestccol)=inhead(n5,htestccol)+hpottydeal(2*n5a,2)-hpottydeal(2*n5a-1,2);
                end
                
                inhead(rachelf+2:nr0,htestccol)=inhead(rachelf+1,htestccol);
                guesserr(nn)=(inhead(rachelf+2,htestccol)-inhead(rachelf+2,hpascol))^2;%error is difference between real and predicted lake hydropotential
                
            end%now we have several steps
            if holyshit==1
                break
            end
            %use newton ralfson method to calcualte next guess
            fprime=(guesserr(3)-guesserr(1))/(2*guessstep);
            f2prime=(guesserr(1)+guesserr(3)-2*guesserr(2))/(guessstep^2);
           % disp([size(hpottydeal),size(rachelsg)]);
            %disp(hpottydeal);
            if abs(f2prime)<=fprime*1e-5%in other word is the 2nd derrivitive flat for some reason?
                Qguessnext=Qguess-0.001;
            else
                Qguessnext=Qguess-fprime/(f2prime);
            end
%             figure(5)
%             plot(hpottydeal(:,1)/km,hpottydeal(:,2:3)/(Rau_w*g));
%             figure(6)
%             plot(inhead(:,scol)/km,[inhead(:,[htestccol,hpasccol,hpascol])/(Rau_w*g),inhead(:,hcol)]);
%             title(staghead(rachelf,Qccol));
%             legend('htest','h-cond','hsheet','h orig');
%             disp([Qguess,Qguessnext,itercount,satis,abs(Qguess-Qguessnext)]);
            
            itercount=itercount+1;
%             disp(itercount);
%             ned=input('iter');
            if itercount>100
                disp([Qguess,Qguessnext])
            end
            if abs(Qguess-Qguessnext)<2e-5;
                satis=4;
                %disp('ooh yeah')
                Qguess=Qguessnext;
                %disp(itercount);
                break
                
               % disp([Qguess,Qguessnext]);
            else
                Qguess=Qguessnext;
                staghead(rachelsg,Qccol)=Qcprev;
                inhead(rachel,Nccol)=Nprev;
                staghead(rachelf,Qccol)=Qguess;
            end
            %satis=input('press 4 to be happy');
            %             if abs(Qguess-staghead(rachelf,Qccol))<0.005
            %                 satis=4;
            %             else
           
        end
      
        
        inhead(rachel,hpasccol)=(inhead(rachel,brkcol)+inhead(rachel,hscol))*Rau_w*g+inhead(rachel,pocol)-inhead(rachel,Nccol);
        %disp([ned,dNcDf,dNcDf*ffstep,Nnew,Nprev,Nprev2])
        
        %         %figure(1)
        %plot(inhead(:,scol)/km,[inhead(:,Tsrcol)*1e10,inhead(:,Qbcol)])
        %eq_2 from kl&ng13  on the main grid
        %disp([Ntot;Ntot2])
        %FORMERLY: staghead(rachelsg,Qccol)
        dhpasc4ds=(inhead(rachel0:rachelf+1,hpasccol)-inhead(rachels0:rachelf,hpasccol))./(inhead(rachel0:rachelf+1,scol)-inhead(rachels0:rachelf,scol));
        % disp(size(dhpasc4ds))
        cutest=sign(dhpasc4ds).*sqrt(abs((dhpasc4ds.*staghead(rachelsg,Sccol).^(8/3))/(fkl*Rau_w*g)));%eq 3 algebraicly reworked
        dQ4now=(staghead(rachel0:rachelf,Qccol)-staghead(rachels0:rachelf-1,Qccol))./(staghead(rachel0:rachelf,scol)-staghead(rachels0:rachelf-1,scol));
        N_alt=0.1*abs(nthroot(max(real(((dQ4now-0*inhead(rachel,mccol)*(1/Rau_w-1/Rau_i)-inhead(rachel,Tsrcol))...
            ./((K0klf)*inhead(rachel,Sccol)))),0),Glen_n));%equation 2 algrbraicly reworked  %note the 0 in front of the term for erositon etc . . 
               
        if stpplot>=1 %flag I had there for the debugging phase of development you will likely never need this, but I keep it here anyway in case something goes horribly wrong
            figure(7)
            plot(inhead(:,scol)/km,[inhead(:,Sccol),inhead_old(:,Sccol),dheaddt(:,Sccol)*1e5,inhead(:,Tsrcol)*1e5],...
                inhead(rachel,scol)/km,[50*inhead(rachel,mccol)/Rau_i,-K0klf*inhead(rachel,Sccol).*inhead(rachel,Nccol).^Glen_n]*1e5,...
                staghead(rachelsg,scol)/km,[50*staghead(rachelsg,mccol)/Rau_i,-K0klf*staghead(rachelsg,Sccol).*staghead(rachelsg,Nccol).^Glen_n]*1e5);
            axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,min([min(inhead(:,Tsrcol))*dtCFLmax*100,max(inhead(:,Sccol))*-.6]),max(inhead(:,Sccol))*1.1]);
            legend('S','S-old','dS','Tsr*100','melt','creep','sgmelt','sgcreep')
            %figure(8)
            %plot(staghead(rachelsg,scol)/km,[staghead(rachelsg,Qccol),Qcprev2]);
            %title('Qc')
            %figure(9)
            %plot(inhead(rachel,scol)/km,[inhead(rachel,Nccol),Nprev2]);
            %title('N')
            figure(8)
            plot(inhead(rachel,scol)/km,[inhead(rachel,Nccol),N_alt/1])
            figure(9)
            plot(staghead(:,scol)/km,-Qstag,inhead(:,scol)/km,inhead(:,Tsrcol)*1e5,...
                staghead(rachelsg,scol)/km,[-staghead(rachelsg,Qccol),-cutest]);
            legend('Qsht','Tsr','QR','Qtest');
            axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,-2.5,max([Qnet,max(staghead(:,Qccol)*1.5)])]);
            figure(10)
            plot(inhead(:,scol)/km,[inhead(:,[hpasccol,hpascol])/(Rau_w*g),inhead(:,hcol)],...
                staghead(rachelsg,scol)/km,staghead(rachelsg,hpasccol)/(Rau_w*g)+staghead(rachelsg,Qccol),...
                staghead(rachelsg,scol)/km,-staghead(rachelsg,Qccol)+min(inhead(:,hcol)),...
                staghead(:,scol)/km,-Qstag+min(inhead(:,hcol)),...
                inhead(rachel,scol)/km,inhead(rachel,Tsrcol)*km+min(inhead(:,hcol)),...
                inhead(:,scol)/km,[min(inhead(:,hcol))+(inhead(:,mccol)*(1/Rau_w-1/Rau_i)+...
                K0kl*inhead(:,Sccol).*inhead(:,Nccol).^Glen_n)*1e5,min(inhead(:,hcol))*ones(nr0,1)]);
            title('hpot')
            legend('h_R','h_s','h_0','HpasC-Q_c','Q_c','Qsht','Tsr','other dQ')
            axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,min(inhead(:,hcol))-5,max(inhead(rachel,hcol))+5]);
            disp(dtCFLmax)
            ned=input('push the button');
            %disp('is this on?')
        end
        %       % break
        %adding material to solve the boundary value problem from the other
        %end
        ih_postiter=inhead;
        sh_postiter=staghead;
        %inhead=ih_preiter;
        %staghead=sh_preiter;
        itercount=0;
    end
    if tsat0+t4now/secyr >= tstop0
        [~,ncVQt]=size(Vl_t);
        %Vl_t(:,trec1+1)=[tsat0+t4now/secyr;inlake(:,Lvolcol)];
        %disp('vl_t6')
    disp(size(Vl_t))
        Qcl_t(:,trec1+1)=[tsat0+t4now/secyr;Qstag(rachelf-1);staghead(rachelf,Qccol)];
        if trec1+2<=ncVQt
            %TSRt6,SSt6,SRt6,hhSt6,hhRt6
            Qcl_t(:,trec1+2:ncVQt)=[];
            Vl_t(:,trec1+2:ncVQt)=[];
            Nt(:,trec1+2:ncVQt)=[];
            TSRt(:,trec1+2:ncVQt)=[];
            SSt(:,trec1+2:ncVQt)=[];
            SRt(:,trec1+2:ncVQt)=[];
            hhSt(:,trec1+2:ncVQt)=[];
            hhRt(:,trec1+2:ncVQt)=[];
        end
        disp('times up donnie')
        break
    end
end
%disp('vl_t7')
%disp(size(Vl_t))
        
disp('years passed')
disp(t4now/secyr);
find(isnan(dQstag)==1)
if plotzon>=1
    figure(5)
    plot(Qcl_t(1,:),Qcl_t(2:3,:));
    overhead52
    figure(6)
    plot(Vl_t(1,:),Vl_t(2:nlake+1,:)/1e9);
    legend('sheetflow','conduitflow');
    overhead52
%      figure(7)
%      plot(staghead(:,scol)/km,[staghead(:,Sccol)],...
%          inhead(:,scol)/km,[inhead_old(:,Sccol),dheaddt(:,Sccol)*5e4,inhead(:,Tsrcol)*5e4],...
%          inhead(rachel,scol)/km,[50*inhead(rachel,mccol)/Rau_i,-K0klf*inhead(rachel,Sccol).*inhead(rachel,Nccol).^Glen_n]*5e4,...
%          staghead(rachelsg,scol)/km,[50*staghead(rachelsg,mccol)/Rau_i,-K0klf*staghead(rachelsg,Sccol).*staghead(rachelsg,Nccol).^Glen_n]*5e4);
%      axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,min([min(inhead(:,Tsrcol))*dtCFLmax*100,max(inhead(:,Sccol))*-.1]),max(inhead(:,Sccol))*1.1]);
%      legend('S','S-old','dS','Tsr*100','melt','creep','sgmelt','sgcreep')
%       figure(8)
%      plot(inhead(rachel,scol)/km,inhead(rachel,mccol)/Rau_s,'b',inhead(rachel,scol)/km,inhead(rachel,mc2col)/Rau_i,'r',...
%          staghead(rachelsg,scol)/km,staghead(rachelsg,mccol)/Rau_s,'b',staghead(rachelsg,scol)/km,staghead(rachelsg,mc2col)/Rau_i,'r',...
%          inhead(rachel,scol)/km,inhead(rachel,Krpccol),'g',staghead(rachelsg,scol)/km,staghead(rachelsg,Krpccol),'g',...
%          inhead(rachel,scol)/km,Krpold,'k',staghead(rachelsg,scol)/km,Krpolds,'k');
%      title('MELT')
%      disp('MELT')
     figure(9)
     plot(inhead(rachel,scol)/km,[inhead(rachel,Vccol),inhead(rachel,Vscol)],staghead(rachelsg,scol)/km,[staghead(rachelsg,Vccol),staghead(rachelsg,Vscol)]);
     title('velocity')
     figure(11)
     plot(inhead(rachel,scol)/km,[inhead(rachel,Cbarccol)],staghead(rachelsg,scol)/km,[staghead(rachelsg,Cbarccol)]);
     title('sed concentration')
     figure(10)
     plot(inhead(:,scol)/km,[inhead(:,[hpasccol,hpascol])/(Rau_w*g),inhead(:,hcol)],...
         staghead(rachelsg,scol)/km,staghead(rachelsg,hpasccol)/(Rau_w*g)+staghead(rachelsg,Qccol),...
         staghead(rachelsg,scol)/km,-staghead(rachelsg,Qccol)+min(inhead(:,hcol)),...
         inhead(rachel,scol)/km,inhead(rachel,Tsrcol)*km+min(inhead(:,hcol)),...
         staghead(:,scol)/km,-Qstag+min(inhead(:,hcol)),...
         inhead(:,scol)/km,[min(inhead(:,hcol))+(inhead(:,mccol)*(1/Rau_w-1/Rau_i)+...
         K0kl*inhead(:,Sccol).*inhead(:,Nccol).^Glen_n)*1e5,min(inhead(:,hcol))*ones(nr0,1)]);
     title('hpot')
     legend('h_R','h_s','h_0','HpasC-Q_c','Q_c','Q_sht','Tsr','other dQ')
     axis([inhead(rachelf,scol)/km-5,inhead(rachel0,scol)/km+5,min(inhead(:,hcol))-5,max(inhead(rachel,hcol))+5]);
%     figure(3)
%     plot(inhead(:,scol)/km,inhead(:,[bedcol,brkcol]),inhead(lkrow0,scol)/km,inhead(lkrow0,[bedcol,brkcol]),'.');
%     figure(4)
%     plot(inhead(:,scol)/km,[inhead(:,hcol),inhead(:,hpascol)/(Rau_w*g),inhead(:,hpascol)/(Rau_w*g)+dheaddt(:,Scol)*1e4,...
%         inhead_preloop(:,hpascol)/(Rau_w*g)]);
%     legend('base','final','+dhdt','step1')
%     figure(5)
%     %had to flip the sign around the darned "s" negative
%     plot(inhead(:,scol)/km,[-inhead(:,[Qbcol]),inhead_preloop(:,Qbcol),-inhead_1(:,Qbcol),-inhead_old(:,Qbcol)],staghead(:,scol)/km,-Qstag);
%     line(inhead(:,scol)/km,-inhead(:,[Qbcol]),'linestyle','none','marker','.')
%     legend('Qb','Qborig','step 1','prevstep','Qstag')
%     figure(6)
%     plot(inhead(:,scol)/km,[inhead(:,dQbcol)])%,dheaddt(:,Scol)]);
%     legend('dQb','dSdt')
%     % ned=input('push a button to end it all now');
end
