function entire
%This program calculates the radius vs. time history of a hole 
%at each of a specifed number of equally spaced depths.  The data
%is also presented as a series of hole profile plots at regular
%intevals before and after the end of drilling.

% The program integrates the partial diff. eq. for
%temperature in the ice.  Each step in the process of
%performing a single time step integration is done on
%a separate array.  Each column in these arrays
%corresponds to a different r location and each row
%starting from the top is a time step.  The order in which
%they are listed is the order in which they are calculated.

%Each section is considered to be independent from those
%below except that the speed at which the lower layers are drilled affects
%the time required to go down to the bottom of the hole and
%then come back to ream.


clc
clear all
delete entiredata.xls
tic %start timer
global rdatprint

format long g
format compact

 

%Starting Parameters 
nlay=25;                    %number of horizontal layers
holedepth=2500;			  %depth of entire hole - meters

%Drill Speeds, in Meters/minute, from bottom layer up
dspd_min(1)=2.25;
dspd_min(2)=2.80;
dspd_min(3)=2.75;
dspd_min(4)=2.60;
dspd_min(5)=2.45;

dspd_min(6)=2.35;
dspd_min(7)=2.20;
dspd_min(8)=2.10;
dspd_min(9)=2.00;
dspd_min(10)=1.90;

dspd_min(11)=1.85;
dspd_min(12)=1.8;
dspd_min(13)=1.75;
dspd_min(14)=1.75;
dspd_min(15)=1.7;

dspd_min(16)=1.7;
dspd_min(17)=1.7;
dspd_min(18)=1.7;
dspd_min(19)=1.7;
dspd_min(20)=1.70;

dspd_min(21)=1.75;
dspd_min(22)=1.75;
dspd_min(23)=1.75;
dspd_min(24)=1.8;
dspd_min(25)=1.8;

%Ream Speeds, in Meters/minute, from bottom layer up
reamspd_min(1)=4.5;
reamspd_min(2)=4.5;
reamspd_min(3)=4.5;
reamspd_min(4)=4.5;
reamspd_min(5)=4.5;

reamspd_min(6)=4.5;
reamspd_min(7)=4.5;
reamspd_min(8)=4.5;
reamspd_min(9)=4.5;
reamspd_min(10)=4.5;

reamspd_min(11)=4.5;
reamspd_min(12)=4.5;
reamspd_min(13)=4.5;
reamspd_min(14)=4.5;
reamspd_min(15)=4.5;

reamspd_min(16)=4.5;
reamspd_min(17)=4.5;
reamspd_min(18)=4.5;
reamspd_min(19)=4.5;
reamspd_min(20)=4.5;

reamspd_min(21)=4.5;
reamspd_min(22)=4.5;
reamspd_min(23)=4.5;
reamspd_min(24)=4.5;
reamspd_min(25)=4.5;
% 


%Spacing of profile plots
befinterval=3600;        %plot interval before end of drilling (sec)
aftinterval=18000;       %plot interval after end of drilling (sec)

%Do Loop Parameters
Tend=270000;

%The following parameters determine the drill tip temp as a function of
%depth
pumprate=1.199e-2;             %flow rate of water down hose - m^3/sec
tswater=80.0;		  %hose temp of water at zero depth (surface) - deg C
%Lambda is the temp drop constant for water travelling in hose
%It is scaled w.r.t its value at 190 gpm (1.199e-2 cu.meters/sec)
lambda=12800*(pumprate/1.199e-2);		  %temp drop time constant for water in hose - n.d.

%Physical Drill Parameters
nozzlediam=0.0254;       %nozzle inside diameter - m
rhose=0.048;                   %outside hose radius - m

%Physical Properties
kice=2.2;                      %thermal conductivity of ice - W/m*C
rhoice=917;                    %density of ice - kg/m^3
%ice density changed from 897 on Jan.6, '05
cice=1950;                     %specific heat of ice - J/kg*C
kwater=0.655;                  %thermal conductivity of water - W/m*C
rhowater=983.3;                  %density of water - kg/m3
%water density changed from 982 on Jan. 6, '05
cwater=4170;                   %specific heat of water - j/kg*C
cfusion=335000;                %feat of fusion of ice - J/kg

%Initial hole radius.  This is the assumed hole radius at the height of
%the tip due to the opening of the hole below the tip.  The resulting hole
%is quite insensitive to the value selected here.  Since the hole is
%opening so fast in this region an error has the effect of translating the
%entire answer up or down a few meters.  A reasonable starting diameter is
%suggested by some papers (Yen and Zehnder, Yen, Gilpin) to be 10 - 15 times
%the nozzle diameter.
rinitial=(10*nozzlediam)/2;         %assumed initial hole radius at tip height - m
%The hose loss coeff. for heat conduction out through the hose 
%is determined by lambda, pumprate, rhowater and cwater - W/m*C
hoseloss=(pumprate*rhowater*cwater)/lambda;		  


%Fill in arrays of variable that change layer to layer

%Calculate the depth of the center of each of the layers and the
%far field temperature there
%(Representative South Pole Ice Temps I've been using:
%-24C at 2250m, -36C at 1750m, -44C at 1250m, -48C at 750m, -50C at
%250m)
laydepth=holedepth/nlay;
for s=1:nlay
    depth(s)=((nlay-s)*laydepth)+(laydepth/2);
    if depth(s)<=1660
        Tinf(s)=(-50.778)+(3*(exp((depth(s)-800)/600)));
    else
        Tinf(s)=(-69.656)+(0.015*depth(s))+(3*(exp((depth(s)-800)/1100)));
    end
end
%Convert dspd_min and reamspd_min to dspd and reamspd (meters/sec)
for s=1:nlay
    dspd(s)=(dspd_min(s))/60;
    reamspd(s)=(reamspd_min(s))/60;
end

%Calculate the times at the center of each layer going
%down and going up.  It is calculated with t=0 when the drill is at the
%bottom of the hole. After drill down time is know the time scale is
%shifted so that t=0 is the time when drilling starts.
downtime=zeros(nlay,1);
uptime=zeros(nlay,1);

downtime(1)=-(laydepth/2)/dspd(1);
uptime(1)=(laydepth/2)/reamspd(1);
for s=2:nlay
    downtime(s)=(downtime(s-1))-(laydepth/2)/dspd(s-1)-(laydepth/2)/dspd(s);
    uptime(s)=(uptime(s-1))+(laydepth/2)/reamspd(s-1)+(laydepth/2)/reamspd(s);
end
Tstartdrill=(downtime(nlay))-(laydepth/2)/dspd(nlay);
Tenddrill=(uptime(nlay))+(laydepth/2)/reamspd(nlay);
%Translate time=0 to start of drilling
for s=1:nlay
    downtime(s)=downtime(s)-Tstartdrill;
    uptime(s)=uptime(s)-Tstartdrill;
end
Tenddrill=Tenddrill-Tstartdrill
Tstartdrill=0;

%Calculated initial values
deltamelt=rhoice/rhowater;

%Set up arrays for collecting profile vs. time data.
nbeforetimes=fix(Tenddrill/befinterval)       %number times before end of drilling
naftertimes=fix((Tend-Tenddrill)/aftinterval)  %number of times after end of drilling
nprofiles=nbeforetimes+naftertimes+1;
%create array of times
lab2='hrs before';
for s=1:(nbeforetimes+1)         
    proftime(s)=Tenddrill+((s-1-nbeforetimes)*befinterval);
    lab1=num2str((Tenddrill-proftime(s))/3600);
    temp=[lab1,lab2];
    proftimehr{s}=(temp);
end

s=nbeforetimes+1;
proftimehr{s}='end of drilling';

lab2='hrs after';
for s=(nbeforetimes+2):nprofiles         
    proftime(s)=Tenddrill+((s-1-nbeforetimes)*aftinterval); 
    lab1=num2str((proftime(s)-Tenddrill)/3600);
    temp=[lab1,lab2];
    proftimehr{s}=(temp);
end
%proftimehr(nbeforetimes+1)='end of drill' hold on this until rest works
% following three lines commented while writing improved section
% profrad=zeros(nlay+1,nprofiles+1);  %initialize data array
% profrad(2:nlay+1,1)=depth ;     %puts depths in first column
% profrad(1,2:nprofiles+1)=proftime;

profrad=cell(nlay+2,nprofiles+1);  %initialize data array
profrad{2,1}='depth(m)';
for n=1:nlay
profrad{n+2,1}=depth(n) ;     %puts depths in first column
end
for n=1:nprofiles
profrad{1,n+1}=proftime(n);
profrad(2,n+1)=proftimehr(n);
end


%Calculate starting water temp in hole
for s=1:nlay
thosewat(s)=tswater*exp(-depth(s)/lambda);		%temp of water inside drill at local depth
%calculate the temp the water has cooled to due to melting ice out to the hole radius
%assumed at drill tip.
tempo1(s)=(rhoice*(cfusion-(cice*Tinf(s))))/(rhowater*cwater);
tempo2(s)=3.14159*(rinitial^2)*dspd(s);
twaterinit(s)=(1/(pumprate+tempo2(s)))*((thosewat(s)*pumprate)-(tempo2(s)*tempo1(s)));
%                				%initial temp of water surrounding tip
%						%after melting out to rinitial - deg. C
%Other calculated values
Qhose(s)=hoseloss*thosewat(s);                     %heat leakage from hose - W/m
Qtip(s)=cwater*rhowater*pumprate*thosewat(s);	%heat flow from tip at local depth - W
end


%Set up the output header and initialize the plot number counter.

%create file that can be indexed through to give start point for plots
plotplaceH=['B2 ';'E2 ';'H2 ';'K2 ';'N2 ';'Q2 ';'T2 ';'W2 ';'Z2 ';'AC2';'AF2';'AI2';...
    'AL2';'AO2';'AR2';'AU2';'AX2';'BA2';'BD2';'BG2';'BJ2';'BM2';'BP2';...
    'BS2';'BV2';'BY2';'CB2';'CE2';'CH2'];
plotplaceV=['B21 ';'E21 ';'H21 ';'K21 ';'N21 ';'Q21 ';'T21 ';'W21 ';'Z21 ';'AC21';'AF21';'AI21';...
    'AL21';'AO21';'AR21';'AU21';'AX21';'BA21';'BD21';'BG21';'BJ21';'BM21';'BP21';...
    'BS21';'BV21';'BY21';'CB21';'CE21';'CH21'];

%create the header that gives info about a plot
headtxt={'Depth','T-inf','T-surf','Dr Spd m/min','Rmspd m/min',...
    'Ream Time','Nozzle diam.','Pump rate','Phase Times','Down f',...
    'Down c','Up','Ream','Freeze','Ream heat per meter',...
    'Ream heat rate'};

%initialize the plot number variable to 1 (first plot)
plotnum=1;




%Ream Loop
reamtime=10000;                 %length of time over which the ream water is assumed to be spread
    
%Step size information    
%radial step size
delr=0.005;                    %initial r step size - m
kdr=1.02;                      %multiplier on r spacing - nd

%Calculate some necessary values
alpha=kice/(rhoice*cice);      %thermal diffusivity
%The following is only constant for constant step ratios in r
rweight=1/(kdr+1);             %the ratio of one r step to two


for s=1:nlay
    %Ream Parameters
    reamspd=reamspd_min(s)/60;			%reamspeed in m/sec
    heatpermeter=Qtip(s)/reamspd;			%ream heat into wall - J/m
    
    %The following definition calculates "reamheat" a value of ream
    %heat that is constant during the ream time.
    reamheat=heatpermeter/reamtime;			%constant heat rate into wall - W/m
    
    %The following is a version that puts the heat into the
    %interface in an exponentially decreasing fashion over the  specified
    %ream time
    ntc=4.0;             %number of time constants over which heat rate decays
    tclength=reamtime/ntc;   %length of one time constant
    HRinit=(heatpermeter)/(tclength*(1-exp(-ntc)));  %initial heat rate beginning of ream
    HR=@(t) HRinit.*exp(-t/tclength);   %heat rate as a function of time - W/m
    %NOTE: The version of the heat rate that gets used is selected in the
    %loop that calls the timestep2 subroutine for reaming.

    Tream=uptime(s)-downtime(s);    %Tream is the time between tip down and tip up



    %Time step size
    delt1=1.0;                     %time step size first part- sec
    if Tream>10000
        nsteps1=10000;                 %number of time steps in first part of drill down
    else
        nsteps1=round(Tream);
    end
    Tstage1=delt1*nsteps1;		    %duration of this part

    delt2=5.0;                                      %time step size second part- sec
    if Tream>10000
        nsteps2=round((Tream-Tstage1)/delt2);                  %number of time steps in second part of drill down
    else
        nsteps2=0;
    end
    Tstage2=delt2*nsteps2;		                     %duration of this part

    %Stage 3 is included if it is desired to treat withdrawing the drill
    %(but before ream) differently from drilling down.  If all time
    %that drillhead is below depth of interest is treated the same then
    %stage 3 is not used (set to zero steps).
    delt3=5.0;                     %time step size third part- sec
    nsteps3=0;                     %number of times steps in third part - from bottom to tip passing going up
    Tstage3=delt3*nsteps3;		    %duration of this part

    delt4=1.0;                     %time step size fourth part (ream)- sec
    nsteps4=round(reamtime/delt4);        %number of times steps in fourth part - ream
    Tstage4=delt4*nsteps4;		    %duration of this part

    delt5=5.0;                                      %time step size fifth part (freeze in) - sec
    nsteps5=round((Tend-downtime(s)-Tstage1-Tstage2-Tstage4)/delt5);   %number of times steps in fifth part - freeze in
    Tstage5=delt5*nsteps5;		                     %duration of this part

    tnst1=nsteps1;
    tnst2=tnst1+nsteps2;
    tnst3=tnst2+nsteps3;
    tnst4=tnst3+nsteps4;
    tnst5=tnst4+nsteps5;
    tnstmax=tnst5;



    %REDIMENSION ARRAYS
    nr=100;                        %number of r values
    nt=tnstmax+2;                    %number of t values
    radiusdat=zeros(500,2);
    r=zeros(nt,nr+1);
    temp=zeros(nt,nr+1);
    dtempdr=zeros(nt,nr+1);
    d2tempdr2=zeros(nt,nr);
    dtempdt=zeros(nt,nr);
    tempeu=zeros(nt,nr+1);
    dtempeu=zeros(nt,nr+1);
    d2tempeu=zeros(nt,nr);
    dtempdteu=zeros(nt,nr);
    tempimpeu=zeros(nt,nr+1);

    clok=zeros(nt+2,1);
    twater=zeros(nt,1);
    diamhyd=zeros(nt,1);
    prandtl=zeros(nt,1);
    visc=zeros(nt,1);
    reynolds=zeros(nt,1);
    hconduct=zeros(nt,1);
    Qinterface=zeros(nt,1);
    Qcond=zeros(nt,1);
    Qmelt=zeros(nt,1);
    dRdt=zeros(nt,1);
    deltar=zeros(nt,1);
    Qholesize=zeros(nt,1);
    Qpastcv=zeros(nt,1);
    iceheat=zeros(nt,1);

    %Fill up initial row of r matrix
    r(1,1) = rinitial;
    varyr=delr;
    for q=2:nr+1
        r(1,q) = (r(1,q-1))+varyr;
        varyr = varyr*kdr;
    end

    %Fill up initial row of temp
    temp(1,1)=0;
    for q=2:nr+1
        temp(1,q)=Tinf(s);
    end

    %set hole water to intitial temperature
    twater(1)=twaterinit(s);


    %set clok to start time
    clok(1)=downtime(s);

    %Do the calculations on qth temps and radii and
    %calculate q+1st temps and radii
    %NOTE: intheat is used only in sub timestep2. It is
    %the heat into the interface.


    %Sect 1) drill, fine time step
    travel=1.2;  %This initializes the travel variable needed to 
    %calculate h in the entrance region.
    delt=delt1;
    for q=1:tnst1
        TIMESTEP1;
    end

    %Sect 2) drill, coarse time step
    delt=delt2;
    for q=(tnst1+1):(tnst2)
        TIMESTEP1;
    end

    %Sect 3) drill being withdrawn, water still flowing (not used)
    delt=delt3;
    intheat=Qhose(s);
    for q=(tnst2+1):(tnst3)
        TIMESTEP2;
    end

    %Sect 4) ream
    delt=delt4;
        %NOTE: intheat must be set to one or the other of the following:
        % If intheat is set to reamheat the reamheat is fed into the
        % interface at a constant average rate over the entire ream time.
        % This is done before entering the loop.
        % If intheat is set to the HR expression the ream heat is fed
        % into the interace at an exponentially decreasing rate over the
        % ream time.  This is done inside the loop since it varies.  
        % In either case the same total heat is fed into the hole wall.
        %intheat=reamheat
    for q=(tnst3+1):(tnst4)
        intheat=0.5*(HR(clok(q-1)-clok(tnst3))+HR(clok(q)-clok(tnst3)));  %This returns the 
        %average value of the heat input at the last ste
        TIMESTEP2;
    end

    %Sect 5) freeze-in
    delt=delt5;
    intheat=0;
    for q=(tnst4+1):(tnst5)
        TIMESTEP2;
        if r(q+1,1)<0.05
            tnst5=q;
            break
        end
    end




    %Print parameters and radius vs. time data

    %Create and print output files to Excel

    headnum=[depth(s),Tinf(s),tswater,dspd_min(s),reamspd_min(s),reamtime,nozzlediam,...
        pumprate,0,Tstage1,Tstage2,Tstage3,Tstage4,Tstage5,heatpermeter,...
        reamheat];
    for q=1:16
        header(q,1)=headtxt(q);
        header(q,2)={headnum(q)};
    end
    header(18,1)={'Time(s)'};
    header(18,2)={'Radius(m)'};
    header(18,3)={'h (W/m^2*C)'};
    %create the column label for this dataset
    lab1=num2str(depth(s));
    lab2='m.deep, ';
    lab3=num2str(dspd_min(s));
    lab4='drill, ';
    lab5=num2str(reamspd_min(s));
    lab6='ream';
    plotname=[lab1,lab2,lab3,lab4,lab5,lab6];
    header(19,2)={plotname};
    xlswrite('entiredata.xls',header,'sheetA',plotplaceH(plotnum,:));
    
    
    
    
    

    %collect the desired data pairs from the clock and radius (clok and r) arrays
    counter=0;
    for q=1:100:5000
        counter=counter+1;
        radiusdat(counter,1)=clok(q);
        radiusdat(counter,2)=r(q,1);
        radiusdat(counter,3)=hconduct(q);
    end

    for q=5001:1000:tnst1
        counter=counter+1;
        radiusdat(counter,1)=clok(q);
        radiusdat(counter,2)=r(q,1);
        radiusdat(counter,3)=hconduct(q);
    end

    for q=tnst1+1:400:tnst2
        counter=counter+1;
        radiusdat(counter,1)=clok(q);
        radiusdat(counter,2)=r(q,1);
        radiusdat(counter,3)=hconduct(q);
    end

    for q=tnst2+1:400:tnst3
        counter=counter+1;
        radiusdat(counter,1)=clok(q);
        radiusdat(counter,2)=r(q,1);
    end

    if tnst1>=5001
        for q=tnst3+1:200:tnst4
            counter=counter+1;
            radiusdat(counter,1)=clok(q);
            radiusdat(counter,2)=r(q,1);
        end
    else
        for q=5001:200:tnst4
            counter=counter+1;
            radiusdat(counter,1)=clok(q);
            radiusdat(counter,2)=r(q,1);
        end
    end

    for q=tnst4+1:400:tnst5+1
        counter=counter+1;
        radiusdat(counter,1)=clok(q);
        radiusdat(counter,2)=r(q,1);
    end

    %establish printing vector of correct length and print to Excel file
    rdatprint=zeros(counter,3);
    %load output data into a file of the correct size
    for q=1:counter
        rdatprint(q,1:3)=radiusdat(q,1:3);
    end
    xlswrite('entiredata.xls',rdatprint,'sheetA',plotplaceV(plotnum,:));
    plotnum=plotnum+1;




    % Collect data for the hole-profile-as-function-of-time plots
    m=1;
    for q=1:nprofiles
        targtime=proftime(q);
        if clok(m)>targtime
            profrad{s+2,q+1}=0;
        else
            while clok(m)<targtime
                m=m+1;
            end
            profrad{s+2,q+1}=r(m,1);
        end
    end



end %end nlay loop


%Print out the column that will plot a r=0.225m line, the
%end-of-drilling line, the 30 hr line and give the end-
%of-drilling time.
targets(1,2)={'radius = 0.225 m'};
targets(2,1)={0};
targets(2,2)={0.225};
targets(3,1)={300000};
targets(3,2)={0.225};

targets(5,2)={'end of drilling'};
targets(6,1)={Tenddrill};
targets(6,2)={0.0};
targets(7,1)={Tenddrill};
targets(7,2)={0.4};

targets(9,2)={'30 hrs after drilling'};
targets(10,1)={Tenddrill+108000};
targets(10,2)={0.0};
targets(11,1)={Tenddrill+108000};
targets(11,2)={0.4};

xlswrite('entiredata.xls',targets,'sheetA',plotplaceH(plotnum,:));
plotnum=plotnum+1;

xlswrite('entiredata.xls',profrad,'sheetA',plotplaceH(plotnum,:));
plotnum=plotnum+1;





% %SUBROUTINES FOLLOW HERE
%
    function TIMESTEP1

        %Calculate qth row of dtempdr
        dtempleft=(temp(q,2)-temp(q,1))/(r(q,2)-r(q,1));
        dtempright=(temp(q,3)-temp(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        dtempdr(q,1)=(dtempleft*(1+rweight))-(dtempright*rweight);
        for k=2:nr
            dtempleft=(temp(q,k)-temp(q,k-1))/(r(q,k)-r(q,k-1));
            dtempright=(temp(q,k+1)-temp(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            dtempdr(q,k)=(dtempleft*(1-rweight))+(dtempright*rweight);
        end
        dtempdr(q,nr+1)=0;

        %Calculate qth row of d2tempdr2
        d2templeft=(dtempdr(q,2)-dtempdr(q,1))/(r(q,2)-r(q,1));
        d2tempright=(dtempdr(q,3)-dtempdr(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        d2tempdr2(q,1)=(d2templeft*(1+rweight))-(d2tempright*rweight);
        for k=2:nr
            d2templeft=(dtempdr(q,k)-dtempdr(q,k-1))/(r(q,k)-r(q,k-1));
            d2tempright=(dtempdr(q,k+1)-dtempdr(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            d2tempdr2(q,k)=(d2templeft*(1-rweight))+(d2tempright*rweight);
        end
        %Calculate the qth row of dtempdt
        for k=1:nr
            dtempdt(q,k)=alpha*(d2tempdr2(q,k)+(dtempdr(q,k)/r(q,k)));
        end

        %Calculate the qth row of tempeu
        tempeu(q,1)=0;
        for k=2:nr
            tempeu(q,k)=temp(q,k)+(delt*dtempdt(q,k));
        end
        tempeu(q,nr+1)=Tinf(s);

        %Calculate qth row of dtempeu (d/dr)
        dtempleft=(tempeu(q,2)-tempeu(q,1))/(r(q,2)-r(q,1));
        dtempright=(tempeu(q,3)-tempeu(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        dtempeu(q,1)=(dtempleft*(1+rweight))-(dtempright*rweight);
        for k=2:nr
            dtempleft=(tempeu(q,k)-tempeu(q,k-1))/(r(q,k)-r(q,k-1));
            dtempright=(tempeu(q,k+1)-tempeu(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            dtempeu(q,k)=(dtempleft*(1-rweight))+(dtempright*rweight);
        end
        dtempeu(q,nr+1)=0;

        %Calculate qth row of d2tempeu (d2/dr2)
        d2templeft=(dtempeu(q,2)-dtempeu(q,1))/(r(q,2)-r(q,1));
        d2tempright=(dtempeu(q,3)-dtempeu(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        d2tempeu(q,1)=(d2templeft*(1+rweight))-(d2tempright*rweight);
        for k=2:nr
            d2templeft=(dtempeu(q,k)-dtempeu(q,k-1))/(r(q,k)-r(q,k-1));
            d2tempright=(dtempeu(q,k+1)-dtempeu(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            d2tempeu(q,k)=(d2templeft*(1-rweight))+(d2tempright*rweight);
        end

        %Calculate the qth row of dtempdteu
        for k=1:nr
            dtempdteu(q,k)=alpha*(d2tempeu(q,k)+(dtempeu(q,k)/r(q,k)));
        end

        %Calculate the qth row of tempimpeu
        tempimpeu(q,1)=0;
        for k=2:nr
            tempimpeu(q,k)=temp(q,k)+(0.5*delt*(dtempdt(q,k)+dtempdteu(q,k)));
        end
        tempimpeu(q,nr+1)=Tinf(s);

        %Calculate change in radius and change in water temperature
        diamhyd(q)=2*(r(q,1)-rhose);
        if twater(q)>=5.8        %find prandtl - approximation blows up towards zero
            prandtl(q)=38*(twater(q)^(-0.6));
        else
            prandtl(q)=13.25;
        end

        if twater(q)>=7.5        %find viscosity - approximation blows up towards zero
            visc(q)=0.006*(twater(q)^(-0.615));
        else
            visc(q)=1.79e-3;
        end

        %reynolds(q)=7.89/(visc(q)*(r(q,1)+rhose));
        %calculate the reynolds number
        flowrate=pumprate+((dspd(s))*pi*rhose^2)-((1-deltamelt)*(dspd(s)*pi*(r(q,1))^2));
        flowveloc=flowrate/(pi*((r(q,1))^2-rhose^2));
        reynolds(q)=(983.3*flowveloc*diamhyd(q))/visc(q);
        %The following selects and entrance region expression for h in the
        %1st section of the calculation at each depth and the selects the
        %fully developed turbulence expression after that.
        if (clok(q)-downtime(s))<20000
            %following is an experimental variation for h calculation based on
            %eqn 6-6 p. 207 of J.P. Holman, Heat Transfer for entrance regions
            hconduct(q)=(kwater/diamhyd(q))*(0.036*(reynolds(q)^0.8)*(prandtl(q)^0.33))*(diamhyd(q)/travel)^0.055;
            travel=travel+(delt*dspd(s));
        else
            %Calculate the boundary layer conduction coefficient from Dittus Boelter
            hconduct(q)=(kwater/diamhyd(q))*(0.023*(reynolds(q)^0.8)*(prandtl(q)^0.3));
        end
             
        Qinterface(q)=2*pi*r(q,1)*hconduct(q)*twater(q);
        Qcond(q)=-2*pi*r(q,1)*kice*0.5*(dtempdr(q,1)+dtempeu(q,1));
        Qmelt(q)=Qinterface(q)-Qcond(q);
        dRdt(q)=Qmelt(q)/(2*pi*r(q,1)*rhoice*cfusion);
        deltar(q)=delt*dRdt(q);
        Qholesize(q)=2*pi*r(q,1)*deltamelt*rhowater*cwater*twater(q)*dRdt(q);
        Qpastcv(q)=(pumprate+(dspd(s)*deltamelt*pi*(r(q,1)^2)))*rhowater*cwater;

        %Calculate radii and temps for the q+1st step
        twater(q+1)=twater(q)+(delt*((Qinterface(q)-Qhose(s)+Qholesize(q))/(-1*Qpastcv(q)/dspd(s))));

        %update clok
        clok(q+1)=clok(q)+delt;

        %Get new hole radius - then adjust radii
        r(q+1,1)=r(q,1)+deltar(q);
        varyr=delr;
        for k=2:nr+1
            r(q+1,k) = (r(q+1,k-1))+varyr;
            varyr = varyr*kdr;
        end

        %Interpolate the temps to match this new grid
        temp(q+1,1)=0;
        if deltar(q)>0
            for k=2:nr
                ratiox=(deltar(q)/(r(q,k+1)-r(q,k)));
                tempx=tempimpeu(q,k+1)-tempimpeu(q,k);
                temp(q+1,k)=tempimpeu(q,k)+ratiox*tempx;
            end
            temp(q+1,nr+1)=Tinf(s);
        elseif deltar(q)<0
            for k=2:nr
                ratiox=(deltar(q)/(r(q,k)-r(q,k-1)));
                tempx=tempimpeu(q,k)-tempimpeu(q,k-1);
                temp(q+1,k)=tempimpeu(q,k)+ratiox*tempx;
            end
            temp(q+1,nr+1)=Tinf(s);
        else                          %for q=0
            for k=1:nr+1
                r(q+1,k)=r(q,k);    %just copying same radii at this time
            end
            temp(q+1,1)=tempimpeu(q,1);
            for k=2:nr
                temp(q+1,k)=tempimpeu(q,k);   %temps the same
            end
            temp(q+1,nr+1)=Tinf(s);
        end
    end %(end of function TIMESTEP1)



    function TIMESTEP2
        % Declare all variables as global
        % global depth Tinf rinitial tswater timeinit reamtime dspd_min
        % global lambda kice rhoice cice kwater rhowater cwater cfusion pumprate rhose hoseloss
        % global thosewat deltamelt drillspeed tempo1 tempo2 twater Qhose Qtip Tend
        % global delr kdr delt1 nsteps1 Tstage1
        % global delt2 nsteps2 Tstage2
        % global delt3 nsteps3 Tstage3
        % global delt4 nsteps4 Tstage4
        % global delt5 nsteps5 Tstage5
        % global tnst1 tnst2 tnst3 tnst4 tnst5 tnstmax
        % global alpha rweight q k
        % global nr nt r temp dtempdr d2tempdr2 dtempdt tempeu dtempeu d2tempeu dtempdteu tempimpeu
        % global clok twater diamhyd Prandtl visc Reynolds hconduct Qinterface Qcond
        % global Qmelt dRdt deltaR Qholesize Qpastcv iceheat varyr delt intheat Tshut
        %Calculate qth row of dtempdr
        dtempleft=(temp(q,2)-temp(q,1))/(r(q,2)-r(q,1));
        dtempright=(temp(q,3)-temp(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        dtempdr(q,1)=(dtempleft*(1+rweight))-(dtempright*rweight);
        for k=2:nr
            dtempleft=(temp(q,k)-temp(q,k-1))/(r(q,k)-r(q,k-1));
            dtempright=(temp(q,k+1)-temp(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            dtempdr(q,k)=(dtempleft*(1-rweight))+(dtempright*rweight);
        end
        dtempdr(q,nr+1)=0;

        %Calculate qth row of d2tempdr2
        d2templeft=(dtempdr(q,2)-dtempdr(q,1))/(r(q,2)-r(q,1));
        d2tempright=(dtempdr(q,3)-dtempdr(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        d2tempdr2(q,1)=(d2templeft*(1+rweight))-(d2tempright*rweight);
        for k=2:nr
            d2templeft=(dtempdr(q,k)-dtempdr(q,k-1))/(r(q,k)-r(q,k-1));
            d2tempright=(dtempdr(q,k+1)-dtempdr(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            d2tempdr2(q,k)=(d2templeft*(1-rweight))+(d2tempright*rweight);
        end

        %Calculate the qth row of dtempdt
        for k=1:nr
            dtempdt(q,k)=alpha*(d2tempdr2(q,k)+(dtempdr(q,k)/r(q,k)));
        end

        %Calculate the qth row of tempeu
        tempeu(q,1)=0;
        for k=2:nr
            tempeu(q,k)=temp(q,k)+(delt*dtempdt(q,k));
        end
        tempeu(q,nr+1)=Tinf(s);

        %Calculate qth row of dtempeu (d/dr)
        dtempleft=(tempeu(q,2)-tempeu(q,1))/(r(q,2)-r(q,1));
        dtempright=(tempeu(q,3)-tempeu(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        dtempeu(q,1)=(dtempleft*(1+rweight))-(dtempright*rweight);
        for k=2:nr
            dtempleft=(tempeu(q,k)-tempeu(q,k-1))/(r(q,k)-r(q,k-1));
            dtempright=(tempeu(q,k+1)-tempeu(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            dtempeu(q,k)=(dtempleft*(1-rweight))+(dtempright*rweight);
        end
        dtempeu(q,nr+1)=0;

        %Calculate qth row of d2tempeu (d2/dr2)
        d2templeft=(dtempeu(q,2)-dtempeu(q,1))/(r(q,2)-r(q,1));
        d2tempright=(dtempeu(q,3)-dtempeu(q,2))/(r(q,3)-r(q,2));
        %rweight=(r(q,2)-r(q,1))/(r(q,3)-r(q,1)) %exclude for contstant step ratio
        d2tempeu(q,1)=(d2templeft*(1+rweight))-(d2tempright*rweight);
        for k=2:nr
            d2templeft=(dtempeu(q,k)-dtempeu(q,k-1))/(r(q,k)-r(q,k-1));
            d2tempright=(dtempeu(q,k+1)-dtempeu(q,k))/(r(q,k+1)-r(q,k));
            %rweight=(r(q,k)-r(q,k-1))/(r(q,k+1)-r(q,k-1)) %exclude for contstant step ratio
            d2tempeu(q,k)=(d2templeft*(1-rweight))+(d2tempright*rweight);
        end

        %Calculate the qth row of dtempdteu
        for k=1:nr
            dtempdteu(q,k)=alpha*(d2tempeu(q,k)+(dtempeu(q,k)/r(q,k)));
        end

        %Calculate the qth row of tempimpeu
        tempimpeu(q,1)=0;
        for k=2:nr
            tempimpeu(q,k)=temp(q,k)+(0.5*delt*(dtempdt(q,k)+dtempdteu(q,k)));
        end
        tempimpeu(q,nr+1)=Tinf(s);


        %CALCULATION OF h DEACTIVATED FOR SHUTOFF AND REAM PROCESSES
        %    %Calculate change in radius and change in water temperature
        %        diamhyd(q)=2*(r(q,1)-rhose)
        %        IF twater(q)>=5.8 then        %find prandtl - approximation blows up towards zero
        %           prandtl(q)=38*(twater(q)^(-0.6))
        %        ELSE
        %           prandtl(q)=13.25
        %        END IF
        %        IF twater(q)>=7.5 then        %find viscosity - approximation blows up towards zero
        %           visc(q)=0.006*(twater(q)^(-0.615))
        %        ELSE
        %           visc(q)=1.79e-3
        %        END IF
        %        reynolds(q)=7.89/(visc(q)*(r(q,1)+rhose))
        %        hconduct(q)=(kwater/diamhyd(q))*(0.023*(reynolds(q)^0.8)*(prandtl(q)^0.3))
        %for reaming reamheat is applied directly to interface
        Qinterface(q)=intheat;
        Qcond(q)=-2*pi*r(q,1)*kice*0.5*(dtempdr(q,1)+dtempeu(q,1));
        Qmelt(q)=Qinterface(q)-Qcond(q);
        dRdt(q)=Qmelt(q)/(2*pi*r(q,1)*rhoice*cfusion);
        deltar(q)=delt*dRdt(q);
        Qholesize(q)=2*pi*r(q,1)*rhowater*cwater*twater(q)*dRdt(q);
        Qpastcv(q)=(pumprate+(dspd(s)*pi*(r(q,1)^2-rhose^2)))*rhowater*cwater;

        %Calculate radii and temps for the q+1st step - held constant for reaming
        %ream heat is added direct to interface not through water temperature
        twater(q+1)=0;   %WATER TEMP HELD CONSTANT DURING SHUTOFF AND REAM
        %This is not accurate as the water must have some nonzero temp in
        %order to pass heat to the ice.  However when the water is still
        %it is hard to evalutate h so the heat generated is passed directly
        %to the interface.  Although the mechanism is inaccurate the result
        %should be realistic.

        %update clok
        clok(q+1)=clok(q)+delt;

        %Get new hole radius - then adjust radii
        r(q+1,1)=r(q,1)+deltar(q);
        varyr=delr;
        for k=2:(nr+1)
            r(q+1,k) = (r(q+1,k-1))+varyr;
            varyr = varyr*kdr;
        end

        %Interpolate the temps to match this new grid
        temp(q+1,1)=0;
        if deltar(q)>0
            for k=2:nr
                ratiox=(deltar(q)/(r(q,k+1)-r(q,k)));
                tempx=tempimpeu(q,k+1)-tempimpeu(q,k);
                temp(q+1,k)=tempimpeu(q,k)+ratiox*tempx;
            end
            temp(q+1,nr+1)=Tinf(s);
        elseif deltar(q)<0
            for k=2:nr
                ratiox=(deltar(q)/(r(q,k)-r(q,k-1)));
                tempx=tempimpeu(q,k)-tempimpeu(q,k-1);
                temp(q+1,k)=tempimpeu(q,k)+ratiox*tempx;
            end
            temp(q+1,nr+1)=Tinf(s);
        else                          %for q=0
            for k=1:nr+1
                r(q+1,k)=r(q,k);    %just copying same radii at this time
            end
            temp(q+1,1)=tempimpeu(q,1);
            for k=2:nr
                temp(q+1,k)=tempimpeu(q,k);   %temps the same
            end
            temp(q+1,nr+1)=Tinf(s);
        end
    end %(end of function TIMESTEP2)
    %
    toc %end timer%
end %end of entire program
