HFBULKTC: computes sensible and latent heat fluxes and other variables.
A=HFBULKTC(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw) computes the following:
Hs = sensible heat flux INTO ocean [W/m^2]
Hl = latent heat flux INTO ocean [W/m^2]
Hl_webb = Webb correction to latent heat flux INTO ocean [W/m^2]
stress = wind stress [N/m^2]
U_star = velocity friction scale [m/s]
T_star = temperature scale [deg C]
Q_star = humidity scale [kg/kg]
L = Monin-Obukhov length [m]
zetu = zr/L
CD = drag coefficient
CT = temperature transfer coefficient (Stanton number)
CQ = moisture transfer coefficient (Dalton number)
RI = bulk Richardson number
Dter = cool-skin temperature difference (optional output) [C];
positive if surface is cooler than bulk (presently no
warm skin permitted by model)
Based on the following buoy input data:
ur = wind speed [m/s] measured at height zr [m]
Ta = air temperature [C] measured at height zt [m]
rh = relative humidity [%] measured at height zq [m]
Pa = air pressure [mb]
Ts = sea surface temperature [C]
sal = salinity [psu (PSS-78)]
(optional - only needed for cool-skin)
dlw = downwelling (INTO water) longwave radiation [W/m^2]
(optional - only needed for cool-skin)
dsw = measured insolation [W/m^2]
(optional - only needed for cool-skin)
nsw = net shortwave radiation INTO the water [W/m^2]
(optional - only needed for cool-skin)
where ur, Ta, rh, Pa, Ts, zr, zt, and zq (and optional sal, dlw,
dsw, and nsw) may be either row or column vectors; and rh, Pa,
zr, zt, and zq (and optional sal) may also be fixed scalars.
Output variables are given as column vectors in A:
1) without cool-skin correction:
A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI]
2) with cool-skin correction:
A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter];0001 function A=hfbulktc(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw) 0002 % HFBULKTC: computes sensible and latent heat fluxes and other variables. 0003 % A=HFBULKTC(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw) computes the following: 0004 % 0005 % Hs = sensible heat flux INTO ocean [W/m^2] 0006 % Hl = latent heat flux INTO ocean [W/m^2] 0007 % Hl_webb = Webb correction to latent heat flux INTO ocean [W/m^2] 0008 % stress = wind stress [N/m^2] 0009 % U_star = velocity friction scale [m/s] 0010 % T_star = temperature scale [deg C] 0011 % Q_star = humidity scale [kg/kg] 0012 % L = Monin-Obukhov length [m] 0013 % zetu = zr/L 0014 % CD = drag coefficient 0015 % CT = temperature transfer coefficient (Stanton number) 0016 % CQ = moisture transfer coefficient (Dalton number) 0017 % RI = bulk Richardson number 0018 % Dter = cool-skin temperature difference (optional output) [C]; 0019 % positive if surface is cooler than bulk (presently no 0020 % warm skin permitted by model) 0021 % 0022 % Based on the following buoy input data: 0023 % 0024 % ur = wind speed [m/s] measured at height zr [m] 0025 % Ta = air temperature [C] measured at height zt [m] 0026 % rh = relative humidity [%] measured at height zq [m] 0027 % Pa = air pressure [mb] 0028 % Ts = sea surface temperature [C] 0029 % sal = salinity [psu (PSS-78)] 0030 % (optional - only needed for cool-skin) 0031 % dlw = downwelling (INTO water) longwave radiation [W/m^2] 0032 % (optional - only needed for cool-skin) 0033 % dsw = measured insolation [W/m^2] 0034 % (optional - only needed for cool-skin) 0035 % nsw = net shortwave radiation INTO the water [W/m^2] 0036 % (optional - only needed for cool-skin) 0037 % 0038 % where ur, Ta, rh, Pa, Ts, zr, zt, and zq (and optional sal, dlw, 0039 % dsw, and nsw) may be either row or column vectors; and rh, Pa, 0040 % zr, zt, and zq (and optional sal) may also be fixed scalars. 0041 % 0042 % Output variables are given as column vectors in A: 0043 % 0044 % 1) without cool-skin correction: 0045 % 0046 % A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI] 0047 % 0048 % 2) with cool-skin correction: 0049 % 0050 % A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter]; 0051 0052 % Code follows Edson and Fairall TOGA COARE code (version 2.0), modified 0053 % to include Rogers' weighting factor for unstable conditions. Code does 0054 % include gustiness, and assumes that the marine boundary layer height is 0055 % known and constant over time for simiplicity. zr/L is limited to 0056 % be <=3.0 to ensure that the code converges to nonzero stress and heat 0057 % flux values for strongly stable conditions. The bulk Richardson number 0058 % is computed between the sea surface and zr as a diagnostic about whether 0059 % turbulent boundary layer theory is applicable. Code does not include 0060 % warm layer effects to modify Ts. See Fairall et al (1996), J. Geophys. 0061 % Res., 101, 3747-3764, for description of full TOGA COARE code and 0062 % comparison with data. 0063 0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0065 % 8/19/98: version 1.1 (rewritten by RP to remove inconsistencies in 0066 % virtual and real temperatures, improve loop structure, 0067 % correct gustiness component of stress computation) 0068 % 4/9/99: version 1.2 (rewritten by AA to clarify some variable names 0069 % and include cool-skin effect and Webb correction to latent 0070 % heat flux added to output matrix) 0071 % 8/5/99: version 2.0 0072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0073 0074 M=length(ur); 0075 0076 % change to column vectors 0077 ur = ur(:); 0078 Ta = Ta(:); 0079 rh = rh(:); 0080 Ts = Ts(:); 0081 Pa = Pa(:); 0082 zr = zr(:); 0083 zt = zt(:); 0084 zq = zq(:); 0085 0086 % create vectors for rh, Pa, zr, zt, and zq, if scalars are input 0087 if length(rh)==1 & M>1 0088 rh=rh*ones(M,1); 0089 end 0090 if length(Pa)==1 & M>1 0091 Pa=Pa*ones(M,1); 0092 end 0093 if length(zr)==1 & M>1 0094 zr=zr*ones(M,1); 0095 end 0096 if length(zt)==1 & M>1 0097 zt=zt*ones(M,1); 0098 end 0099 if length(zq)==1 & M>1 0100 zq=zq*ones(M,1); 0101 end 0102 0103 if nargin > 8 0104 % optional cool-skin stuff 0105 sal = sal(:); 0106 dlw = dlw(:); 0107 dsw = dsw(:); 0108 nsw = nsw(:); 0109 % create vector for sal if scalar is input 0110 if length(sal)==1 & M>1 0111 sal=sal*ones(M,1); 0112 end 0113 end 0114 0115 % initialize various constants 0116 as_consts; 0117 0118 tol=.001; % tolerance on Re changes to make sure soln has converged. 0119 0120 onethird=1./3; 0121 o61=1/eps_air-1; % 0.61 (moisture correction for temperature) 0122 0123 visc=viscair(Ta); % viscosity 0124 Qsats=Qsat_coeff*qsat(Ts,Pa); % saturation specific humidity; the Qsat_coeff 0125 % value is set in routine as_consts.m 0126 Q=(0.01.*rh).*qsat(Ta,Pa); % specific humidity of air [kg/kg] 0127 T =Ta+CtoK; % convert to K 0128 Tv=T.*(1 + o61*Q); % air virtual temperature 0129 rho=(100*Pa)./(gas_const_R*Tv); % air density 0130 Dt=(Ta+0.0098.*zt)-Ts; % adiabatic temperature difference 0131 Dq=Q-Qsats; % humidity difference 0132 0133 % compute initial neutral scaling coefficients 0134 S=sqrt(ur.^2 + min_gustiness.^2); 0135 cdnhf=sqrt(cdntc(S,zr,Ta)); % Smith's neutral cd as first guess 0136 0137 z0t=7.5*10^(-5); 0138 ctnhf=kappa./log(zt./z0t); 0139 0140 z0q=z0t; 0141 cqnhf=kappa./log(zq./z0q); 0142 0143 U_star = cdnhf.*S; % (includes gustiness) 0144 T_star = ctnhf.*Dt; % 0145 Q_star = cqnhf.*Dq; % 0146 0147 Dter = 0; 0148 Dqer = 0; 0149 if nargin > 8 0150 % initial cool-skin thickness guess 0151 delta = 0.001*ones(size(Ts)); 0152 end 0153 0154 Reu=0;Ret=0;Req=0; 0155 % begin iteration loop to compute best U_star, T_star, and Q_star 0156 for iter1=1:80; 0157 0158 ReuO=Reu; RetO=Ret; ReqO=Req; % Save old values 0159 0160 % Compute Monin-Obukov length (NB - definition given as eqn (7) 0161 % of Fairall et al (1996) probably wrong, following, e.g. 0162 % Godfrey and Bellars (1991), JGR, 96, 22043-22048 and original code) 0163 bs=g*(T_star.*(1 + o61*Q) + o61*T.*Q_star)./Tv; 0164 L=(U_star.^2)./(kappa*bs); 0165 % set upper limit on zr/L = 3.0 to force convergence under 0166 % very stable conditions. Assume that zr, zt and zq comparable. 0167 index_limit = (L<zr/3 & L>0); 0168 L(index_limit)=zr(index_limit)/3; 0169 0170 zetu=zr./L; % nondimensionalized heights 0171 zett=zt./L; 0172 zetq=zq./L; 0173 0174 % surface roughness 0175 z0=(Charnock_alpha/g).*U_star.^2 + R_roughness.*visc./U_star; 0176 0177 % compute U_star correction for non-neutral conditions 0178 cdnhf=kappa./(log(zr./z0)-psiutc(zetu)); 0179 U_star=cdnhf.*S; 0180 0181 Reu=z0.*U_star./visc; % roughness Reynolds # 0182 [Ret,Req]=LKB(Reu); % compute other roughness Reynolds #s 0183 0184 % compute t and q roughness scales from roughness R#s 0185 z0t=visc.*Ret./U_star; 0186 z0q=visc.*Req./U_star; 0187 0188 % compute new transfer coefficients at measurement heights 0189 cthf=kappa./(log(zt./z0t)-psittc(zett)); 0190 cqhf=kappa./(log(zq./z0q)-psittc(zetq)); 0191 0192 % compute new values of T_star, Q_star 0193 T_star=cthf.*(Dt + Dter); 0194 Q_star=cqhf.*(Dq + Dqer); 0195 0196 % estimate new gustiness 0197 Ws=U_star.*(-CVB_depth./(kappa*L)).^onethird; 0198 wg=min_gustiness*ones(M,1); 0199 j=find(zetu<0); % convection in unstable conditions only 0200 wg(j)=max(min_gustiness,beta_conv.*Ws(j)); % set minimum gustiness 0201 S=sqrt(ur.^2 + wg.^2); 0202 0203 if nargin > 8 0204 % compute cool-skin parameters 0205 [delta,Dter,Dqer] = cool_skin(sal,Ts-Dter,rho,cp,Pa, ... 0206 U_star,T_star,Q_star, ... 0207 dlw,dsw,nsw,delta,g,gas_const_R, ... 0208 CtoK,Qsat_coeff); 0209 end 0210 0211 end % end of iteration loop 0212 0213 ii= abs(Reu-ReuO)>tol | abs(Ret-RetO)>tol | abs(Req-ReqO)>tol; 0214 if any(ii), 0215 disp(['Algorithm did not converge for ' int2str(sum(ii)) ' values. Indices are:']); 0216 disp(find(ii)'); 0217 warning('Not converged!'); 0218 end; 0219 0220 0221 % compute latent heat 0222 Le=(2.501-0.00237*(Ts-Dter))*10^6; 0223 0224 % compute fluxes into ocean 0225 Hs=rho.*cp.*U_star.*T_star; 0226 Hl=rho.*Le.*U_star.*Q_star; 0227 0228 % compute transfer coefficients at measurement heights 0229 CD=(U_star./S).^2; 0230 CT=U_star.*T_star./(S.*(Dt + Dter)); % Stanton number 0231 CQ=U_star.*Q_star./(S.*(Dq + Dqer)); % Dalton number 0232 0233 % to compute mean stress, we don't want to include the effects 0234 % of gustiness which average out (in a vector sense). 0235 stress=rho.*CD.*S.*ur; 0236 0237 % compute bulk Richardson number (as a diagnostic) - the "T" 0238 % is probably not quite right - assumes T \ approx. Ts (good enough though) 0239 RI=g.*zr.*((Dt + Dter) + o61*T.*(Dq + Dqer))./(Tv.*S.^2); 0240 0241 % compute Webb correction to latent heat flux into ocean 0242 W = 1.61.*U_star.*Q_star + (1 + 1.61.*Q).*U_star.*T_star./T; % eqn. 21 0243 Hl_webb = rho.*Le.*W.*Q; % eqn. 22, Fairall et al. (1996), JGR, 101, p3751. 0244 0245 % output array 0246 if nargin > 8 0247 % output additional cool-skin parameters 0248 A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter]; 0249 else 0250 % otherwise 0251 A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI]; 0252 end 0253 0254 0255 0256 function y=psiutc(zet) 0257 % PSIUTC: computes velocity profile function following TOGA/COARE. 0258 % y=PSIUTC(zet) computes the turbulent velocity profile function given 0259 % zet = (z/L), L the Monin-Obukoff length scale, following Edson and 0260 % Fairall TOGA COARE code (version 2.0) as modified to include Rogers' 0261 % weighting factor to combine the Dyer and free convection forms for 0262 % unstable conditions. 0263 0264 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0265 % 8/28/98: version 1.1 0266 % 8/5/99: version 1.2 0267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0268 0269 c13=1.0./3.0; 0270 sq3=sqrt(3.0); 0271 0272 % stable conditions 0273 y=-4.7*zet; 0274 0275 % unstable conditions 0276 j=find(zet<0); 0277 zneg=zet(j); 0278 0279 % nearly stable (standard functions) 0280 x=(1-16.0.*zneg).^0.25; 0281 y1=2.0.*log((1+x)./2) + log((1+x.^2)./2) -2.*atan(x) + pi/2; 0282 0283 % free convective limit 0284 x=(1-12.87*zneg).^c13; 0285 y2=1.5*log((x.^2+x+1)./3) - sq3*atan((2.*x+1)/sq3) + pi/sq3; 0286 0287 % weighted sum of the two 0288 F=1.0./(1+zneg.^2); 0289 y(j)=F.*y1+(1-F).*y2; 0290 0291 0292 0293 0294 function y=psittc(zet) 0295 % PSITTC: computes potential temperature profile following TOGA/COARE. 0296 % y=PSITTC(zet) computes the turbulent potential temperature profile 0297 % function given zet = (z/L), L the Monin-Obukoff length scale, following 0298 % Edson and Fairall TOGA COARE code (version 2.0), as modified to use 0299 % Rogers' weighting factor to combine the Dyer and free convective 0300 % forms for unstable conditions. 0301 0302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0303 % 8/28/98: version 1.1 0304 % 8/5/99: version 1.2 0305 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0306 0307 c13=1.0./3.0; 0308 sq3=sqrt(3.0); 0309 0310 % stable conditions 0311 y=-4.7.*zet; 0312 0313 % unstable conditions 0314 j=find(zet<0); 0315 zneg=zet(j); 0316 0317 % nearly stable (standard functions) 0318 x=(1-16.0.*zneg).^0.25; 0319 y1=2.0*log((1+x.^2)./2); 0320 0321 % free convective limit 0322 x=(1-12.87*zneg).^c13; 0323 y2=1.5.*log((x.^2+x+1)./3.0) - sq3.*atan((2.*x+1)./sq3) + pi./sq3; 0324 0325 % weighted sum of the two 0326 F=1.0./(1+zneg.^2); 0327 y(j)=F.*y1 + (1-F).*y2; 0328 0329 0330 0331 0332 0333 function [Ret,Req]=LKB(Reu); 0334 % LKB: computes rougness Reynolds numbers for temperature and humidity 0335 % [Ret,Req]=LKB(Reu) computes the roughness Reynolds for temperature 0336 % and humidity following Liu, Katsaros and Businger (1979), J. Atmos. 0337 % Sci., 36, 1722-1735. 0338 0339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0340 % 8/28/98: version 1.1 0341 % 8/5/99: version 1.2 0342 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0343 0344 Ret=.177*ones(size(Reu)); 0345 Req=.292*ones(size(Reu)); 0346 j=find(Reu>.11 & Reu<=.825); 0347 Ret(j)=1.376.*Reu(j).^0.929; 0348 Req(j)=1.808.*Reu(j).^0.826; 0349 j=find(Reu>.825 & Reu<=3); 0350 Ret(j)=1.026./Reu(j).^0.599; 0351 Req(j)=1.393./Reu(j).^0.528; 0352 j=find(Reu>3 & Reu<=10); 0353 Ret(j)=1.625./Reu(j).^1.018; 0354 Req(j)=1.956./Reu(j).^0.870; 0355 j=find(Reu>10 & Reu<=30); 0356 Ret(j)=4.661./Reu(j).^1.475; 0357 Req(j)=4.994./Reu(j).^1.297; 0358 j=find(Reu>30); 0359 Ret(j)=34.904./Reu(j).^2.067; 0360 Req(j)=30.790./Reu(j).^1.845;