% Script for 1st order closure model for CBL with no
% shear.

%  Input parameters and theta profile (thetai):
%  Use evenly spaced flux levels from 0-ztop m with spacing dz, with 
%  initial theta given at half-levels zh. lambda is the asymptotic Blackadar 
%  lengthscale, and Hs0 is the surface sensible heat flux (W/m2)

  tf = 3*3600;
  dtout = 600;
  dz = 100;
  ztop = 1000;   
  lambda = 100;
  Hs0 = 300;
  zh = (dz/2:dz:(ztop-dz/2))';
  thetai = 290 + 0.01*zh;

  tspan = 0:dtout:tf;
  zf = (0:dz:ztop)';
  nz = length(zh);

  rhoref = 1.2;
  Cp = 1000;
  Fth0 = Hs0/(rhoref*Cp);

%  Specify Blackadar master lengthscale for interior flux points for use in
%  ODE solver.

  k = 0.4;
  l = lambda./(1 + lambda./(k*zf(2:nz)));

% Use Matlab ODE solver to solve the equation. Stiff solver works better!

  [t,theta] = ode15s('dthdt',tspan,thetai,[],zf,zh,l,Fth0);
  theta = theta';  % Now each column of theta corresponds to a time in tspan.

  subplot(2,2,1)
  colormap(gray)
  pcolor(t/3600,zh,theta)
  shading interp
  hold on
  contour(t/3600,zh,theta,10,'k')
  colorbar
  xlabel('time (hr)')
  ylabel('z (m)')
  axis([0 tf/3600 0 ztop])
  title('Time-height section of \theta')
  hold off

  subplot(2,2,2)
  ihrs = (mod(t,3600)-1)<1;  % Pick out indices for hourly values
  plot(theta(:,ihrs),zh)
  xlabel('\theta (K)')
  ylabel('z (m)')
  title('Hourly profiles of \theta')
  legend('0 hr','1 hr','2 hr','3 hr')
  text(0.5 + min(thetai), 0.75*ztop, ['\lambda = ' num2str(lambda) ' m'])
  text(0.5 + min(thetai), 0.9*ztop, ['H_s(0) = ' num2str(Hs0) ' W m^{-2}'])

%  Generate hourly profiles of heat flux and K

  for hr = 1:3
    thetahr = theta(:,ihrs);    
    [tend, K, Fth] = dthdt(hr*3600,thetahr(:,hr+1),flag,zf,zh,l,Fth0);
    Khr(:,hr+1) = [0; K; 0];
    Fthhr(:,hr+1) = Fth;
  end
  subplot(2,2,3)
  plot(Fthhr,zf)
  xlabel('Sensible heat flux (W m^{-2})')
  ylabel('z (m)')
  title('Hourly heat flux profiles')
  legend('0 hr','1 hr','2 hr','3 hr')

  subplot(2,2,4)
  plot(Khr,zf)
  xlabel('K (m s^{-1})')
  ylabel('z (m)')
  title('Hourly K profiles')
  legend('0 hr','1 hr','2 hr','3 hr')



  
