function [gamma_th, b1, b2, b3, b4] = cbl_coeffs(mon, ws, h, um0, beta)
% function to calculate the ambient stability of the cbl model
% as well as the coefficients b1 to b4, 
% from the month, wind speed, cbl height, mixed-layer u and beta

if (h <= 50)
 switch mon
  case {6,7,8}			%winter
   if ws <= 5  gamma_th = 200e-3;
   else gamma_th = 25e-3;
   end;
  case {9,10,11}		%spring
   if ws <= 5  gamma_th = 50e-3;
   else gamma_th = 10e-3;
   end;
  case {12,1,2}			%summer
   if ws <= 5  gamma_th = 10e-3;
   else gamma_th = 6.7e-3;
   end;
  case {3,4,5}			%autumn
   if ws <= 5  gamma_th = 150e-3;
   else gamma_th = 25e-3;
   end;
 end;

elseif (h>50 & h<=100)
  switch mon
  case {6,7,8}			%winter
   if ws <= 5  gamma_th = 9.2e-3;
   else gamma_th = 25e-3;
   end;
  case {9,10,11}		%spring
   if ws <= 5  gamma_th = 25e-3;
   else gamma_th = 10e-3;
   end;
  case {12,1,2}			%summer
   if ws <= 5  gamma_th = 10e-3;
   else gamma_th = 6.7e-3;
   end;
  case {3,4,5}			%autumn
   if ws <= 5  gamma_th = 9.2e-3;
   else gamma_th = 25e-3;
   end;
 end;

elseif (h > 100)
  switch mon
  case {6,7,8}			%winter
   gamma_th = 9.2e-3;
  case {9,10,11}		%spring
   gamma_th = 9.2e-3;
  case {12,1,2}			%summer
   gamma_th = 6.7e-3;
  case {3,4,5}			%autumn
   gamma_th = 9.2e-3;
 end;

end; 	%if statement

%calculate b coefficients
 cp = 1004; 		%specific heat capacity
 b1 = 2*(1+2*beta)./(um0*gamma_th*cp);
 b2 = (1+beta)./(um0*cp);
 b3 = um0*gamma_th*cp./(2*(1+2*beta));
 b4 = (um0*cp)./(1+beta);