% Supplemental program 6.1% -------------------------------------------------------------------------% Calculate friction velocity and sensible heat flux given wind speed and% temperature at two heights using Monin-Obukhov similarity theory or% roughness sublayer theory from Physick and Garratt (1995)% -------------------------------------------------------------------------% --- Physical constantsrgas=8.31446;% Universal gas constant (J/K/mol)var.k=0.4;% von Karman constantvar.g=9.80665;% Gravitational acceleration (m/s2)cpair=29.2;% Specific heat of air at constant pressure (J/mol/K)% --- Input variablesvar.d=19.0;% Displacement height (m)var.z1=21.0;% Height (m)var.u1=1.0;% Wind speed at height z1 (m/s)var.t1=29.0+273.15;% Temperature at height z1 (K)var.z2=29.0;% Height (m)var.u2=2.1;% Wind speed at height z2 (m/s)var.t2=28.1+273.15;% Temperature at height z2 (K)var.zstar=49.0;% Height of roughness sublayer (m)abl='MOST';% Use Monin-Obukhov similarity theory% abl = 'RSL'; % Use roughness sublayer theory% --- Molar density (mol/m3)rhomol=101325/(rgas*var.t2);switchabl% -----------------------------------% Monin-Obukhov similarity theory% -----------------------------------case'MOST'% Use bisection to solve for L as specified by the function "most"% and then calculate fluxes for that value of LL1=100;% Initial guess for Obukhov length L (m)L2=-100;% Initial guess for Obukhov length L (m)func_name='most';% The function name is "most", in the file most.m[L]=bisect(func_name,L1,L2,0.01,var);% Solve for L (m)% Evaluate psi for momentum and scalars at heights z2 and z1[psi_m_z2]=psi_m_monin_obukhov((var.z2-var.d)/L);[psi_m_z1]=psi_m_monin_obukhov((var.z1-var.d)/L);[psi_c_z2]=psi_c_monin_obukhov((var.z2-var.d)/L);[psi_c_z1]=psi_c_monin_obukhov((var.z1-var.d)/L);% Calculate u* and T* and the sensible heat fluxustar=(var.u2-var.u1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1));tstar=(var.t2-var.t1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1));H=-rhomol*cpair*tstar*ustar;% Calculate aerodynamic conductancesgam=rhomol*var.k*ustar/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1));gac=rhomol*var.k*ustar/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1));fprintf('Monin-Obuhkov similarity theory\n')fprintf('L = %15.3f\n',L)fprintf('u* = %15.3f\n',ustar)fprintf('T* = %15.3f\n',tstar)fprintf('H = %15.3f\n',H)fprintf('gam = %15.3f\n',gam)fprintf('gac = %15.3f\n',gac)% -----------------------------------% Roughness sublayer theory% -----------------------------------case'RSL'% Use bisection to solve for L as specified by the function "rsl"% and then calculate fluxes for that value of LL1=100;% Initial guess for Obukhov length L (m)L2=-100;% Initial guess for Obukhov length L (m)func_name='rsl';% The function name is "rsl", in the file rsl.m[L]=bisect(func_name,L1,L2,0.01,var);% Solve for L (m)% Evaluate psi for momentum and scalars at heights z2 and z1[psi_m_z2]=psi_m_monin_obukhov((var.z2-var.d)/L);[psi_m_z1]=psi_m_monin_obukhov((var.z1-var.d)/L);[psi_c_z2]=psi_c_monin_obukhov((var.z2-var.d)/L);[psi_c_z1]=psi_c_monin_obukhov((var.z1-var.d)/L);% Evaluate the roughness sublayer-modified psi (between z1 and z2)f1_psi_m_rsl=@(z)(1-16*(z-var.d)/L).^(-0.25).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f1_psi_c_rsl=@(z)(1-16*(z-var.d)/L).^(-0.50).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f2_psi_m_rsl=@(z)(1+5*(z-var.d)/L).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f2_psi_c_rsl=@(z)(1+5*(z-var.d)/L).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);if(L<0)psi_m_rsl=integral(f1_psi_m_rsl,var.z1,var.z2);psi_c_rsl=integral(f1_psi_c_rsl,var.z1,var.z2);elsepsi_m_rsl=integral(f2_psi_m_rsl,var.z1,var.z2);psi_c_rsl=integral(f2_psi_c_rsl,var.z1,var.z2);end% Calculate u* and T* and the sensible heat fluxustar=(var.u2-var.u1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1)-psi_m_rsl);tstar=(var.t2-var.t1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1)-psi_c_rsl);H=-rhomol*cpair*tstar*ustar;% Calculate aerodynamic conductancesgam=rhomol*var.k*ustar/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1)-psi_m_rsl);gac=rhomol*var.k*ustar/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1)-psi_c_rsl);fprintf('Roughness sublayer theory\n')fprintf('L = %15.3f\n',L)fprintf('u* = %15.3f\n',ustar)fprintf('T* = %15.3f\n',tstar)fprintf('H = %15.3f\n',H)fprintf('gam = %15.3f\n',gam)fprintf('gac = %15.3f\n',gac)end
function[c]=bisect(func_name,a,b,delta,var)% -----------------------------------------------------------------% Use the bisection method to find the root of a function f% between a and b. The root is refined until its accuracy is delta.%% Input: func_name ! Name of the function to solve% a ! Low endpoint of the interval% b ! High endpoint of the interval% delta ! Tolerance/accuracy% var ! Input variables for function% Output: c ! Root% -----------------------------------------------------------------% Evaluate function at a and bfa=feval(func_name,a,var);fb=feval(func_name,b,var);% Error check: root must be bracketedif(sign(fa)==sign(fb))error('bisect error: f must have different signs at the endpoints a and b')end% Iterate to find rootwhile(abs(b-a)>2*delta)c=(b+a)/2;fc=feval(func_name,c,var);if(sign(fc)~=sign(fb))a=c;fa=fc;elseb=c;fb=fc;endend
function[fx]=most(x,var)% -------------------------------------------------------------------------% Use Monin-Obukhov similarity theory to obtain the Obukhov length (L).%% This is the function to solve for the Obukhov length. For current estimate% of the Obukhov length (x), calculate u* and T* and then the new length (L).% The function value is the change in Obukhov length: fx = x - L.%% Input: x ! Current estimate for Obukhov length (m)% var.z1 ! Height (m)% var.z2 ! Height (m)% var.u1 ! Wind speed at z1 (m/s)% var.u2 ! Wind speed at z2 (m/s)% var.t1 ! Temperature at z1 (m/s)% var.t2 ! Temperature at z2 (m/s)% var.d ! Displacement height (m)% var.k ! von Karman constant% var.g ! Gravitational acceleration (m/s2)% Output: fx ! Change in Obukhov length (x - L)%% Local: psi_m_z2 ! psi for momentum at height z2 (dimensionless)% psi_m_z1 ! psi for momentum at height z1 (dimensionless)% psi_c_z2 ! psi for scalars at height z2 (dimensionless)% psi_c_z1 ! psi for scalars at height z1 (dimensionless)% ustar ! Friction velocity (m/s)% tstar ! Temperature scale (K)% L ! Obukhov length (m)% -------------------------------------------------------------------------% Prevent near-zero values of Obukhov lengthif(abs(x)<=0.1)x=0.1;end% Evaluate psi for momentum at heights z2 and z1[psi_m_z2]=psi_m_monin_obukhov((var.z2-var.d)/x);[psi_m_z1]=psi_m_monin_obukhov((var.z1-var.d)/x);% Evaluate psi for scalars at heights z2 and z1[psi_c_z2]=psi_c_monin_obukhov((var.z2-var.d)/x);[psi_c_z1]=psi_c_monin_obukhov((var.z1-var.d)/x);% Calculate u* (m/s) and T* (K)ustar=(var.u2-var.u1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1));tstar=(var.t2-var.t1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1));% Calculate L (m)L=ustar^2*var.t2/(var.k*var.g*tstar);% Calculate change in Lfx=x-L;
function[psi_c]=psi_c_monin_obukhov(x)% --- Evaluate the Monin-Obukhov psi function for scalars at xif(x<0)y=(1-16*x)^0.25;psi_c=2*log((1+y^2)/2);elsepsi_c=-5*x;end
function[psi_m]=psi_m_monin_obukhov(x)% --- Evaluate the Monin-Obukhov psi function for momentum at xif(x<0)y=(1-16*x)^0.25;psi_m=2*log((1+y)/2)+log((1+y^2)/2)-2*atan(y)+pi/2;elsepsi_m=-5*x;end
function[fx]=rsl(x,var)% -------------------------------------------------------------------------% Use Physick and Garratt (1995) roughness sublayer theory (RSL) to% obtain the Obukhov length (L).%% This is the function to solve for the Obukhov length. For current estimate% of the Obukhov length (x), calculate u* and T* and then the new length (L).% The function value is the change in Obukhov length: fx = x - L.%% Input: x ! Current estimate for Obukhov length (m)% var.z1 ! Height (m)% var.z2 ! Height (m)% var.u1 ! Wind speed at z1 (m/s)% var.u2 ! Wind speed at z2 (m/s)% var.t1 ! Temperature at z1 (m/s)% var.t2 ! Temperature at z2 (m/s)% var.d ! Displacement height (m)% var.k ! von Karman constant% var.g ! Gravitational acceleration (m/s2)% var.zstar ! Height of roughness sublayer (m)% Output: fx ! Change in Obukhov length (x - L)%% Local: psi_m_z2 ! psi for momentum at height z2 (dimensionless)% psi_m_z1 ! psi for momentum at height z1 (dimensionless)% psi_c_z2 ! psi for scalars at height z2 (dimensionless)% psi_c_z1 ! psi for scalars at height z1 (dimensionless)% psi_m_rsl ! roughness sublayer-modified psi for momentum (dimensionless)% psi_c_rsl ! roughness sublayer-modified psi for scalars (dimensionless)% ustar ! Friction velocity (m/s)% tstar ! Temperature scale (K)% L ! Obukhov length (m)% -------------------------------------------------------------------------% Prevent near-zero values of Obukhov lengthif(abs(x)<=0.1)x=0.1;end% Evaluate psi for momentum at heights z2 and z1[psi_m_z2]=psi_m_monin_obukhov((var.z2-var.d)/x);[psi_m_z1]=psi_m_monin_obukhov((var.z1-var.d)/x);% Evaluate psi for scalars at heights z2 and z1[psi_c_z2]=psi_c_monin_obukhov((var.z2-var.d)/x);[psi_c_z1]=psi_c_monin_obukhov((var.z1-var.d)/x);% Evaluate the roughness sublayer-modified psi (between z1 and z2)f1_psi_m_rsl=@(z)(1-16*(z-var.d)/x).^(-0.25).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f1_psi_c_rsl=@(z)(1-16*(z-var.d)/x).^(-0.50).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f2_psi_m_rsl=@(z)(1+5*(z-var.d)/x).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);f2_psi_c_rsl=@(z)(1+5*(z-var.d)/x).*(1-exp(-0.7*(1-(z-var.d)/(var.zstar-var.d))))./(z-var.d);if(x<0)psi_m_rsl=integral(f1_psi_m_rsl,var.z1,var.z2);psi_c_rsl=integral(f1_psi_c_rsl,var.z1,var.z2);elsepsi_m_rsl=integral(f2_psi_m_rsl,var.z1,var.z2);psi_c_rsl=integral(f2_psi_c_rsl,var.z1,var.z2);end% Calculate u* (m/s) and T* (K)ustar=(var.u2-var.u1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_m_z2-psi_m_z1)-psi_m_rsl);tstar=(var.t2-var.t1)*var.k/(log((var.z2-var.d)/(var.z1-var.d))-(psi_c_z2-psi_c_z1)-psi_c_rsl);% Calculate L (m)L=ustar^2*var.t2/(var.k*var.g*tstar);% Calculate change in Lfx=x-L;