% cbl_growth.m
% a model of convective boundary layer growth, following Garratt Ch6, 
% with adaptions by Renfrew and King, Feb 99
% 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 (Smith 88, sfcflux_all.m)
% hsh(x) = sensible heat flux from height eqn.
% hst(x) = sensible heat flux from theta eqn.
clear all

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

%-----------------------------------------------------------------------
% input data - at distance x=0. Ropex = 1, Chang&Braham = 2, Labsea = 3
% test case = 4, mizclm (Guest et al.) = 5, mefis (heinemann) = 6
%-----------------------------------------------------------------------
input_set = 6

if (input_set == 1) 
% Ropex cruise data from 4 Feb 98 - input data from mpaws
% variables are ropex_thetam0, ropex_pm0, ropex_um0; times 1:24 hrs 4th Feb
 load U:/data/matlab/mpaws_cbl_input;
 xmax = 40;			%fetch max (km)
 beta = [0.2] 			%entrainment parameter
 gamma = 6.4e-3; 		%ambient stability K/m - from Ropex radiosondes
 h0 = [100]			%cbl height (m)
 thetam0 = [ropex_thetam0(11)]	%mixed layer theta from AWS (K)
 ts0 = [273];			%sea surface temp (from cruise) (K) 
 pm0 = [ropex_pm0(11)]		%mixed layer pressure from AWS (mb)
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [ropex_um0(11)]		%mixed layer wind speed from AWS (m/s)
 z = [3 h0 2 40];		%measurement heights for bulk flux routine

elseif (input_set ==2)
% Chang and Braham (1991, MWR) Lake Michigan validation data
 xmax = 150;			%fetch max (km)
 beta = [0.2] 			%entrainment parameter
 gamma = 10.0e-3; 		%ambient stability K/m 
 h0 = [200]			%cbl height 
 thetam0 = [246.8];		%chang and braham (fig5)
 ts0 = [274];			%sst chang and braham = 1C
 pm0 = [995];			%changed as heights are 50m asl
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
% u30 = [9.14];			%mean of winds in table 4
 u30 = [10];			%mean quoted on page 2272
 z = [50 50 50 0];		%measurement heights for bulk flux routine
	
elseif (input_set ==3)
% Renfrew and Moore (1999) Labsea validation data from 8th Feb
 xmax = 400;			%fetch max (km)
 beta = [0.1] 			%entrainment parameter
 gamma = 5.57e-3; 		%ambient stability K/m - labsea estimate
 h0 = [200]			%cbl height 
 thetam0 = [254.317];		%labsea
 ts0 = [274];			
 pm0 = [1000];			
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [13.449];
 z = [100 100 100 0];		%measurement heights for bulk fluxes

elseif (input_set ==4)
% test case, used for paper
 xmax = 50;			%fetch max (km)
 beta = [0.2] 			%entrainment parameter
% gamma = 3.9e-3; 		%ambient stability K/m - summer min 
% 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 = [50]			%cbl height 
 thetam0 = [253];		
 ts0 = [274];			
 pm0 = [1000];			
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [5];
 z = [3 50 2 40];		%measurement heights for bulk fluxes

elseif (input_set ==5)
% mizclm data set of spring arctic radiosonde data from Guest et al. 
 xmax = 150;			%fetch max (km)
 beta = [0.2]; 			%entrainment parameter
 gamma = 5.0e-3; 		%ambient stability K/m - estimated
 h0 = [475];			%cbl height 
 thetam0 = [259.5];		
 ts0 = [272];			
 pm0 = [1000];			
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [7.55];
 z = [10 10 10 0];		%measurement heights for bulk fluxes

elseif (input_set ==6)
% mefis data set from Filchner/Ronne Heinemann (1988) paper & tethersondes
 xmax = 50;			%fetch max (km)
 beta = [0.2]; 			%entrainment parameter
 gamma = 11.4e-3; 		%ambient stability K/m 
 h0 = [175];			%cbl height 
 thetam0 = [260.5];		
 ts0 = [271];			
 pm0 = [991];			
 [thetas0, thetaes0, qs0, qsats0, qsatis0] = thermo_rh(ts0-273.15,pm0,100);
 u30 = [2];
 z = [50 50 50 0];		%measurement heights for bulk fluxes
end;

% 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
%old heat flux calculation, with ctheta constant
% ctheta = 1.4e-3; 		%heat coefficient constant
% hs0 = rhom0*cp*ctheta*u30*(thetas0-thetam0); 	%heat flux (W/m2)
% [hs0, lhf0] = sfcflux_hf(u30, tm0-273, rhm0, ts0-273, pm0, z); 

% use robust surface flux routine------------
% [hs0,lhf0,hf,wstar,tau,cd,ce,cdn,ch,ustar,tstar,qstar,uh0,L]=...
% sfcflux_all(u30,tm0-273,rhm0,ts0-273,pm0,z,h0); %fluxes using u30 at ht z(1)
 [hs0,lhf0,hf,wstar,tau,cd,ce,cdn,ch,ustar,tstar,qstar,uh0,L]=...
 sfcflux_all_robust(u30,tm0-273,rhm0,ts0-273,pm0,z,h0); %fluxes using u30 at ht z(1)

%mixed layer wind speed, calcuated from sfcflux routine for height h0
 um0 = uh0;	

%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
% [hsb(1), lhf(1)] = sfcflux_hf(u30, tm(1)-273, rhm0, ts0-273, pm0, z); 
 [hsb(1), lhf(1)] = sfcflux_hf_robust(u30, tm(1)-273, rhm0, ts0-273, pm0, z); 

%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));	
% [hsb(1), lhf(1)] = sfcflux_hf(u30, tm(1)-273, rhm0, ts0-273, pm0, z); 
 [hsb(1), lhf(1)] = sfcflux_hf_robust(u30, tm(1)-273, rhm0, ts0-273, pm0, z); 
  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)); 	
%  [hsb(i), lhf(i)] = sfcflux_hf(u30, tm(i)-273, rhm0, ts0-273, pm0, z); 
  [hsb(i), lhf(i)] = sfcflux_hf_robust(u30, tm(i)-273, rhm0, ts0-273, pm0, z); 
   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 U:/data/matlab/cbl_growth  x h thetam hs 



