% cbl_growth_gm.m
% a model of convective boundary layer growth, following Garratt Ch6, 
% with adaptions by Renfrew and King, Feb 99
% simplified for the UEA Geophysical Modelling Course.
% 
% Assumes steady state solution and estimates cbl growth with fetch, x.
% h(x) = cbl height, 
% thetam(x) = mixed layer potential temp, 
% hs(x)  = sensible heat flux from bulk formula (simplified to constant C_H)
% hsh(x) = sensible heat flux from height eqn.
% hst(x) = sensible heat flux from theta eqn.
clear all

%constants
R = 287;		        %specifc gas constant
cp = 1004; 	            %specific heat capacity
K= R/cp;

%-----------------------------------------------------------------------
% test case - see Renfrew and King 2000, figure 3
% set upstream (ice shelf) input variables and parameters
 xmax = 50;			    %fetch max (km)
 beta = [0.2] 			%entrainment parameter
% gamma = 6.4e-3; 		%ambient stability K/m - winter min
 gamma = 9.2e-3; 		%ambient stability K/m - winter mean 
% gamma = 12.0e-3; 		%ambient stability K/m - winter max
 h0 = [100]			    %cbl height 
 thetam0 = [253];		
 ts0 = [271];			
 pm0 = [1000];	
 % calculate thermodynamic variables, assume RH w.r.t. ice is 100%
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [5];

% common input variables
a = ((pm0./1000).^K);
tm0 = thetam0.*a;		    %mixed layer temperature (K)
rhm0 = 80;			        %mixed layer rh (%) 
rhom0 = (pm0*100)./(R*tm0);	%mixed layer rho using surface p

% heat flux calculation, with ctheta constant
 ctheta = 1.4e-3; 		                        % heat coefficient constant
 hs0 = rhom0*cp*ctheta*u30*(thetas0-thetam0); 	% heat flux (W/m2)

%mixed layer wind speed, set as u30 in this case
um0 = u30;	

%---------------------------------
%integration set up
%---------------------------------
deltax= 0.5			            %constant deltax (km) (default = 0.1)
x = [deltax:deltax:xmax]*1e3; 	%distance (m)
hepsilon = 0.01; 		        %height convergence must be to hepsilon (m)
hseps = 2;			            %heat flux convergence check (W/m2)

%iteration coefficients - note um0 used for model growth
b1 = 2*(1+2*beta)/(um0*gamma*cp);
b2 = (1+beta)/(um0*cp);
b3 = um0*gamma*cp/(2*(1+2*beta));
b4 = (um0*cp)/(1+beta);

%--------------------------------------------
%initial step, ie. i=1
%--------------------------------------------
i=1;
dx(1) = x(1)-0;
h(1) = sqrt(h0*h0 + b1*dx(1)*(hs0/rhom0));
thetam(1) = thetam0 + b2*dx(1)*(hs0/(rhom0*h0));
tm(1) = thetam(1)*a; 
rhom(1) = (pm0*100)/(R*tm(1));	
hsh(1) = b3*rhom(1)*(h(1)*h(1) - h0*h0)/dx(1);		%hs from height eqn
hst(1) = b4*rhom(1)*h(1)*(thetam(1)-thetam0)/dx(1);	%hs from theta eqn
hsb(1) = rhom(1)*cp*ctheta*u30*(thetas0-thetam(1)); 	%hs from bulk formula

%--------------------------------------------
%iteration loop for i=1
%--------------------------------------------
oldh= h0;
iter =0;
hs(1) = hsb(1);		%set hs to hsb		
while (abs(h(1)-oldh) > hepsilon)
    oldh = h(1);
    h(1) = sqrt(h0*h0 + b1*(dx(1)/2)*(hs0/rhom0 + hs(1)/rhom(1)));
    thetam(1) = thetam0 + b2*(dx(1)/2)*(hs0/(rhom0*h0) + hs(1)/(rhom(1)*h(1)));
    tm(1) = thetam(1)*a; 
    rhom(1) = (pm0*100)/(R*tm(1));	
    hsh(1) = b3*rhom(1)*(h(1)*h(1) - h0*h0)/dx(1);
    hst(1) = b4*rhom(1)*h(1)*(thetam(1)-thetam0)/dx(1);
    hsb(1) = rhom(1)*cp*ctheta*u30*(thetas0-thetam(1));	
    hs(1) = hsb(1);		%set hs to hsb		
    iter = iter + 1;
end;

%check hs are sufficiently close
if ((abs(hsb(1)-hsh(1)) > hseps) | (abs(hsb(1)-hst(1)) > hseps))
    disp('warning hs values differ by more than'); hseps
end;

%--------------------------------------------
%integration over fetch ie. i = 2:length(x) 
%including an iteration for convergence of equations
%--------------------------------------------
for i = 2:1:length(x),
    dx(i) = x(i)-x(i-1);
    h(i) = h(i-1); rhom(i) = rhom(i-1); thetam(i) = thetam(i-1); 	%first guess 
    hs(i) = hs(i-1);		%set hs to hs from previous x position
    oldh = h0;
    iter = 0;
    while (abs(h(i)-oldh) > hepsilon)
        oldh = h(i);
        h(i) = sqrt(h0*h0 + b1*(dx(i)/2)* ...
            (hs0/rhom0 + hs(i)/rhom(i) + 2*sum(hs(1:i-1)./rhom(1:i-1))));
        thetam(i) = thetam0 + b2*(dx(i)/2)* ...
            (hs0/(rhom0*h0) + hs(i)/(rhom(i)*h(i)) + ...
            2*sum(hs(1:i-1)./(rhom(1:i-1).*h(1:i-1))));
        tm(i) = thetam(i)*a; 
        rhom(i) = (pm0*100)/(R*tm(i));	
        hsh(i) = b3*rhom(i)*(h(i)*h(i) - h(i-1)*h(i-1))/dx(i);
        hst(i) = b4*rhom(i)*h(i)*(thetam(i)-thetam(i-1))/dx(i);
        hsb(i) = rhom(i)*cp*ctheta*u30*(thetas0-thetam(i)); 	
        hs(i) = hsb(i); 		
        iter = iter + 1;
    end;
    
    %check hs are sufficiently close
    if ((abs(hsb(i)-hsh(i)) > hseps) | (abs(hsb(i)-hst(i)) > hseps))
        disp('warning hs values differ by more than');
        hseps
        i
    end;
end; %for x loop

%add in shore point
x = [0 x./1000]; 		%fetch in km
h = [h0 h];    		%height (m)
thetam = [thetam0 thetam];	%theta (K)
hs = [hs0 hs];			%sensible heat flux (W/m2)
rhom = [rhom0 rhom];		%density
tm = [tm0 tm];			%temperature (K)

save C:/data/uea/env3a04/cbl_growth  x h thetam hs