z2=10;z1=2;
u2=5;u1=4;t2=0;t1=-1;
th=t1+273.15;
ri=9.82/th*(t2-t1)/(u2-u1).^2*(z2-z1);%1-level bulk gradie
zm=sqrt(z2*z1);%geometric mean height
zms=(z2-z1)/log(z2/z1);%analytical neutral height

L0=1*zm/ri;
L=L0;
%%%%%%%  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
[ri L sqrt(z2*z1)/L L*ri/sqrt(z1*z2)]

lz=-1:.01:3;
z=exp(lz);
zr=zm;
    fm=phiu_36(z/L);%phi function as derivative of psi function
    ft=phit_36(z/L);
    Ri=z/L.*ft./fm.^2;%MO similarity ri as function of z/L
ii=find(abs(Ri)>abs(ri));%find value of z corrresponding too Rig=Ri
j=ii(1)-1;   
    
[ri Ri(j) z(j) zm zms 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;title(['z2= ' num2str(z2) ' z1= ' num2str(z1) ' u2= ' num2str(u2) ' u1= ' num2str(u1) ' th2= ' num2str(t2) ' th1= ' num2str(t1) ' Ri= ' num2str(ri) ' 1/L= ' num2str(1/L)]) 


us=.4*(u2-u1)/(log(z2/z1)-psiu_30(z2/L)+psiu_30(z1/L));
ts=.4*(t2-t1)/(log(z2/z1)-psit_30(z2/L)+psit_30(z1/L));
u10n=u2+us/.4*(log(10/z2)+psiu_30(z2/L));

[us ts .4*z(j)*9.83/th*ts/us^2 1e3*(us/u10n)^2 -1.3*1e3*us*ts]

