Home > air_sea > hfbulktc.m

hfbulktc

PURPOSE ^

HFBULKTC: computes sensible and latent heat fluxes and other variables.

SYNOPSIS ^

function A=hfbulktc(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw)

DESCRIPTION ^

 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];

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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;

Generated on Thu 26-Apr-2007 12:23:13 by m2html © 2003