% Compute bulk Richardson number and adjusted height for Ri-bulk
% 12/2021 CF

z2=10; z1=2;    % level heights, m
u2=5; u1=4;     % wind speeds
t2=0; t1=-1;    % temperature
th = t1+273.15;
ri = 9.82/th*(t2-t1)/(u2-u1).^2*(z2-z1);  % 1-level bulk gradient Ri
zm = sqrt(z2*z1);         % geometric mean estimated height of Ri
zms = (z2-z1)/log(z2/z1); % analytical neutral height of Ri

L = 1*zm/ri;  % first guess for L, using zm
%%%%%%%  iterate to find value of L consistent with ri
for i=1:50
    Lm = ri/(z2-z1)*(log(z2/z1)-psiu_36(z2/L)+psiu_36(z1/L))^2/(log(z2/z1)-psit_36(z2/L)+psit_36(z1/L));
    L = 1/Lm;
end
display(['Ri = ',sprintf('%6.3f',ri),';  L = ',sprintf('%6.3f',L),...
         ';  zm/L = ',sprintf('%6.3f',sqrt(z2*z1)/L),...
         ';  L*ri/sqrt(z1*z2) = ',sprintf('%6.3f',L*ri/sqrt(z1*z2))])

% find value of z corrresponding to Ri(MOS) = Ri
lz = -1:.01:3;
z = exp(lz);                % z grid, 0.37 to 20m
zr = zm;
[fm,ft] = phiut_36(z/L);
Ri = z/L.*ft./fm.^2;        % MO similarity Ri as function of z/L
ii = find(abs(Ri)>abs(ri));
j = ii(1)-1;                % z(j) = best value of z for computed Ri

display(['Ri = ',sprintf('%6.3f',ri),';  Ri(MOS) = ',sprintf('%6.3f',Ri(j)),...
         ';  z(best) = ',sprintf('%6.3f',z(j)),';  z(geo mean) = ',sprintf('%6.3f',zm),...
         ';  z(neutral) = ',sprintf('%6.3f',zms),';  z(j)/L = ',sprintf('%6.3f',z(j)/L)])

figure;
semilogy(Ri,z,'-',[ri ri],[z(1) z(length(z))],'-',[Ri(1) Ri(length(z))],[z1 z1],'-',...
        [Ri(1) Ri(length(z))],[z2 z2],'-',ri,z(j),'o');
xlabel('Ri(z)'); ylabel('z (m)'); grid on; grid minor;
legend('Ri(z)','Ri_b','z1','z2','Ri(z-best)','Location','southeast');
title(['z2= ' num2str(z2),'  z1= ' num2str(z1),'  u2= ' num2str(u2),...
       '  u1= ' num2str(u1),'  t2= ' num2str(t2),'  t1= ' num2str(t1),...
       '  Ri= ' num2str(ri),'  1/L= ' num2str(1/L)]);

usr = 0.4*(u2-u1)/(log(z2/z1) - psiu_36(z2/L) + psiu_36(z1/L));
tsr = 0.4*(t2-t1)/(log(z2/z1) - psit_36(z2/L) + psit_36(z1/L));
u10n = u2 + usr/0.4*(log(10/z2) + psiu_36(z2/L));

display(['usr = ',sprintf('%6.3f',usr),';  tsr = ',sprintf('%6.3f',tsr),...
         ';  z/L = ',sprintf('%6.3f',0.4*z(j)*9.83/th*tsr/usr^2),...
         ';  Cd = ',sprintf('%6.3f',1e3*(usr/u10n)^2),...
         ';  Hs = ',sprintf('%6.3f',-1.3*1e3*usr*tsr)])
