%% MOex.m:  Monin-Obukhov example script
%  Find surface buoyancy and momentum fluxes given u and thetav
%  at a known reference height, above a surface of known roughness
%  and thetav.

%  User parameters

z0 = 1e-3;   % Surface roughness length [m]
zT0 = 0.15*z0; % Typical ratio of thermal to mom. roughness length
thv0 = 300;  % Surface virtual temperature thetav0 [K]
zR = 2;  % Reference height [m]
uR = 5;  % u(zR) [m/s]
thvR = 290;  % thetav(zR) [K]
zmin = 0.1;  % Min height on profile plots
zmax = 30; % Max height on profile plots

%  Physical parameters
k = 0.4; % von Karman constant
g = 9.8; % gravity [m/s2]
beta = 5; % Constant in stable M-O stability functions
gamma = 16; % Constant in unstable M-O stability functions
%%  Iterate between Cd, Ch and zetaR = zR/L to convergence

% Initial guess is a neutral surface layer (zetaR=0,buoy flux B0=0)

zetaR = 0

% Iterate 5 times (hopefully enough to achieve convergence)
% Check convergence by printing out zetaR after each iteration.
for n=1:5
  if(zetaR>=0) % Stable/neutral
    Psim = -beta*zetaR;
    Psih = Psim;
  else  % Unstable
    x = (1 - gamma*zetaR)^0.25;
    fact1 = 0.5*(1+x^2);
    fact2 = 0.25*(1+x)^2;
    Psim = log(fact1*fact2) - 2*atan(x)+pi/2;
    Psih = 2*log(fact1);
  end
  Cd = k^2/(log(zR/z0)-Psim)^2;
  Ch = k^2/((log(zR/z0)-Psim)*(log(zR/zT0)-Psih));
  upwp0 = -Cd*uR^2;
  wpthvp0 = Ch*uR*(thv0-thvR);
  ustar = sqrt(-upwp0);
  B0 = g*wpthvp0/thv0;
  Lob = -ustar^3/(k*B0);
  zetaR = zR/Lob  % Updated zetaR
end

%% Plot wind and theta profiles

log10z = log10(zmin):0.1:log10(zmax);
z = 10.^log10z;
zeta = z/Lob;

if(Lob>=0) % Stable/neutral
  Psim = -beta*zeta;
  Psih = Psim;
else  % Unstable
  x = (1 - gamma*zeta).^0.25;
  fact1 = 0.5*(1+x.^2);
  fact2 = 0.25*(1+x).^2;
  Psim = log(fact1.*fact2) - 2*atan(x)+pi/2;
  Psih = 2*log(fact1);
end

% Wind and theta profiles
u = (ustar/k)*(log(z/z0) - Psim);

thvstar = wpthvp0/ustar;
thv = thv0 - (thvstar/k)*(log(z/zT0) - Psih);

subplot(1,2,1)
semilogy(u,z)
xlabel('u')
ylabel('z [m]')
ylim([zmin zmax])
hold on
plot([u(1) u(end)],[zR zR],'--')
hold off

subplot(1,2,2)
semilogy(thv,z)
xlabel('\theta_v [K]')
ylabel('z [m]')
ylim([zmin zmax])
hold on
plot([thv(1) thv(end)],[zR zR],'--')
hold off
