%%  Script for AS547 HW 3

%% Problem 1

%  Parts (a) and (b)

  z     = [  2       4       8      16      32   ];
  u     = [  5.81    6.70    7.49    8.14    8.66];
  theta = [307.20  306.65  306.28  305.88  305.62];

  g = 9.8;
  theta0 = mean(theta);
  lnz = log(z);
  lnzm = 0.5*(lnz(1:4) + lnz(2:5));
  zm = exp(lnzm)
  dudz = diff(u)./(zm.*diff(lnz))
  dthetadz = diff(theta)./(zm.*diff(lnz))
  dbdz = (g/theta0)*dthetadz;
  Ri = dbdz./(dudz).^2
  L = mean(zm./Ri)
  subplot(2,2,1)
  plot(Ri,zm,'+',[0 zm/L],[0 zm],'-')
  xlabel('Ri')
  ylabel('z_m')

% Part (c): calculate ustar, sensible heat flux using zm(1) = 2.8 m gradients.

  k = 0.4;
  rho = 1.2;
  cp = 1000;
  zeta = zm(1)/L;
  phim = (1 - 16*zeta)^(-1/4);
  ustar = k*zm(1)*dudz(1)/phim
  surfbuoyflux = -ustar^3/(k*L);
  surfheatflux = rho*cp*surfbuoyflux*(theta0/g)  


% Calculate roughness length from z(1) = 2 m data

  x = (1 - 16*z(1)/L)^(1/4);
  psim = log((1+x^2)/2*((1+x)/2)^2) - 2*atan(x) + pi/2;
  lnzz0 = k*u(1)/ustar + psim;
  z0 = z(1)*exp(-lnzz0)

%% Problem 2:  CDN vs. U10 over the oceans from Charnock and bulk formula

k = 0.4;  % von Karman constant
ac = 0.016;  % Charnock coefficient
g = 9.8;  % gravity [m/s2]
u10 = 4:20;   % 10 m wind
z10 = 10;   % 10 m

% Iterative CDN from Charnock

% First guess ustar2 and z0 from bulk CDN

CDNbulk = (0.75+0.067*u10)*1e-3; 

% Now iterate CDN and z0 to convergence
% niter = 5 iterations gets the relative difference between
% successive iterations of CDN below 10^-4 at all wind speeds

CDN = CDNbulk;  % Initial guess at CDN
for i = 1:5
  ustar2 = CDN.*u10.^2;
  z0 = (ac/g)*ustar2;
  CDNold = CDN;
  CDN = k.^2./(log(z10./z0)).^2;
  maxrelerr = max(abs(CDNold./CDN - 1))
end

subplot(2,2,1)
plot(u10,CDN,'b',u10,CDNbulk,'r')
xlabel('u_{10} [m/s]')
ylabel('C_{DN}')
legend('Charnock','Bulk','Location','NorthWest')

%% Problem 3

%  Constants and parameters of problem

  k = 0.4; % von Karman constant
  g = 9.8; % gravity [m/s2]
  cp = 1004; % Isobaric heat capacity of air [J/kg-K]
  Lv = 2.5e6; % Latent heat of vaporization [J/kg]
  SST = 284;  % Sea-surface temperature [K]
  rho = 1.2; % Air density
  z10 = 10;  % Reference height [m]
  u10 = 10; % 10 m wind speed [m/s]
  theta10 = 287; % 10 m potential temperature 
                 % (referenced to surface pressure)

% (a) Neutral friction velocity and surface stress

  Cdn = (0.75 + 0.067*u10)*1e-3
  ustar = sqrt(Cdn*u10^2)
  stress = rho*ustar^2

% (b) Bulk-formula heat fluxes and stability parameter

  q10 = 0.0074;  % 10 m mixing ratio
  qs0 = 0.0079; % Surface saturation mixing ratio
  Chn = 1.3e-3;
  LHF = rho*Lv*Chn*u10*(qs0-q10) % [W/m2]
  SHF = rho*cp*Chn*u10*(SST-theta10) % [W/m2]
  B0 =  g/(rho*cp*theta10) * (SHF + (0.61*cp*theta10/Lv)*LHF)   %[m2/s3]

  Lob = -ustar.^3/(k*B0)
  zetaR = z10/Lob
