diff --git a/README.md b/README.md new file mode 100644 index 0000000..2869677 --- /dev/null +++ b/README.md @@ -0,0 +1 @@ +# air-sea diff --git a/air_dens.m b/air_dens.m new file mode 100644 index 0000000..51c453b --- /dev/null +++ b/air_dens.m @@ -0,0 +1,29 @@ +function rhoa=air_dens(Ta,RH,Pa) +% AIR_DENS: computes computes the density of moist air. +% rhoa=AIR_DENS(Ta,RH,Pa) computes the density of moist air. +% Air pressure is optional. +% +% INPUT: Ta - air temperature Ta [C] +% RH - relative humidity [%] +% Pa - air pressure (optional) [mb] +% +% OUTPUT: rhoa - air density [kg/m^3] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 4/7/99: version 1.2 (contributed by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% load constants +as_consts; + +if nargin == 2 + Pa = P_default; +end + +o61 = 1/eps_air-1; % 0.61 (moisture correction for temp.) +Q = (0.01.*RH).*qsat(Ta,Pa); % specific humidity of air [kg/kg] +T = Ta+CtoK; % convert to K +Tv = T.*(1 + o61*Q); % air virtual temperature +rhoa = (100*Pa)./(gas_const_R*Tv); % air density [kg/m^3] + diff --git a/albedo.m b/albedo.m new file mode 100644 index 0000000..15d6bcc --- /dev/null +++ b/albedo.m @@ -0,0 +1,33 @@ +function alb=albedo(trans,sunalt) +% ALBEDO: computes sea surface albedo following Payne (1972). +% alb=ALBEDO(trans,sunalt) computes the sea surface albedo from the +% atmospheric transmittance and sun altitude by linear interpolation +% using Table 1 in Payne (1972), J. Atm. Sci., 29, 959-970. Assumes +% trans and sunalt both matrices of same size. Table 1 is called +% albedot1.mat. +% +% INPUT: trans - atmospheric transmittance +% sunalt - sun altitude [deg] +% +% OUTPUT: alb - albedo + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/10/96: version 1.0 +% 7/24/98: version 1.1 (rev. to handle out-of-range input values by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% load table 1 +load albedot1 + +% create axes +x=[0:2:90]; +y=[0:.05:1.0]'; + +alb=ones(size(trans))+NaN; +k=sunalt>0 & finite(trans) & trans<=1.01; + +% interpolate +alb(k)=interp2(x,y,albedot1,sunalt(k),trans(k)); + + diff --git a/albedot1.mat b/albedot1.mat new file mode 100644 index 0000000..add83d4 Binary files /dev/null and b/albedot1.mat differ diff --git a/as_consts.m b/as_consts.m new file mode 100644 index 0000000..7af3ce5 --- /dev/null +++ b/as_consts.m @@ -0,0 +1,77 @@ +% AS_CONSTS: returns values of many constants used in AIR_SEA TOOLBOX +% AS_CONSTS: returns values of many constants used in the AIR_SEA +% TOOLBOX. At end of this file are values of constants from COARE +% to be used when running the test program t_hfbulktc.m + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/19/98: version 1.1 (contributed by RP) +% 4/7/99: version 1.2 (revised to include COARE test values by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------- physical constants + +g = 9.8; % acceleration due to gravity [m/s^2] +sigmaSB = 5.6697e-8; % Stefan-Boltzmann constant [W/m^2/K^4] +eps_air = 0.62197; % molecular weight ratio (water/air) +gas_const_R = 287.04; % gas constant for dry air [J/kg/K] +CtoK = 273.16; % conversion factor for [C] to [K] + + +% ------- meteorological constants + +kappa = 0.4; % von Karman's constant +Charnock_alpha = 0.011; % Charnock constant (for determining roughness length + % at sea given friction velocity), used in Smith + % formulas for drag coefficient and also in Fairall + % and Edson. use alpha=0.011 for open-ocean and + % alpha=0.018 for fetch-limited (coastal) regions. +R_roughness = 0.11; % limiting roughness Reynolds # for aerodynamically + % smooth flow + + +% ------ defaults suitable for boundary-layer studies + +cp = 1004.7; % heat capacity of air [J/kg/K] +rho_air = 1.22; % air density (when required as constant) [kg/m^2] +Ta_default = 10; % default air temperature [C] +P_default = 1020; % default air pressure for Kinneret [mbars] +psych_default = 'screen'; % default psychmometer type (see relhumid.m) +Qsat_coeff = 0.98; % satur. specific humidity coefficient reduced + % by 2% over salt water + + +% the following are useful in hfbulktc.m +% (and are the default values used in Fairall et al, 1996) + +CVB_depth = 600; % depth of convective boundary layer in atmosphere [m] +min_gustiness = 0.5; % min. "gustiness" (i.e., unresolved fluctuations) [m/s] + % should keep this strictly >0, otherwise bad stuff + % might happen (divide by zero errors) +beta_conv = 1.25;% scaling constant for gustiness + + +% ------ short-wave flux calculations + +Solar_const = 1368.0; % the solar constant [W/m^2] represents a + % mean of satellite measurements made over the + % last sunspot cycle (1979-1995) taken from + % Coffey et al (1995), Earth System Monitor, 6, 6-10. + + +% ------ long-wave flux calculations + +emiss_lw = 0.985; % long-wave emissivity of ocean from Dickey et al + % (1994), J. Atmos. Oceanic Tech., 11, 1057-1076. + +bulkf_default = 'berliand'; % default bulk formula when downward long-wave + % measurements are not made. + + +% ------ constants used for COARE; to use simply delete the % +% g = 9.7803; % acceleration due to gravity [m/s^2] +% sigmaSB = 5.67e-8; % Stefan-Boltzmann constant [m^2/K^4] +% gas_const_R = 287.1; % gas constant for dry air [J/kg/K] +% cp = 1004.67; % heat capacity of air [J/kg/K] +% beta_conv = 1.20; % scaling constant for gustiness +% emiss_lw = 0.97; % long-wave emissivity diff --git a/binave.m b/binave.m new file mode 100644 index 0000000..447191c --- /dev/null +++ b/binave.m @@ -0,0 +1,50 @@ +function bindat = binave(data,r) +% BINAVE: averages vector data in bins of length r. +% bindat=BINAVE(data,r) computes an average vector of the vector +% data in bins of length r. The last bin may be the average of +% less than r elements. Useful for computing daily average time +% series (with r=24 for hourly data). +% +% INPUT: data - data vector +% r - number of elements in bin to be averaged +% +% OUTPUT: bindat - bin-averaged vector + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 9/19/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% check input +if nargin < 2 + error('Not enough input arguments.') +end +if abs(r-fix(r)) > eps + error('Bin size R must be a positive integer.') +end +if fix(r) == 1 + bindat = data; + return +end +if r <= 0 + error('Bin size R must be a positive integer.') +end + +[N,M]=size(data); + +% compute bin averaged series +l = length(data)/r; +l = fix(l); +bindat = mean(reshape(data(1:l*r),r,l)); + +if length(data)>l*r, + bindat=[bindat,mean(data(l*r+1:end))]; +end; + +if N~=1 + bindat=bindat'; +end + + + diff --git a/blwhf.m b/blwhf.m new file mode 100644 index 0000000..d14039f --- /dev/null +++ b/blwhf.m @@ -0,0 +1,84 @@ +function lwest = blwhf(Ts,Ta,rh,Fc,bulk_f) +% BLWHF: estimates net long-wave heat flux using bulk formulas +% lwest = BLWHF(Ts,Ta,rh,Fc,bulk_f) estimates the net longwave +% heat flux into the ocean using one of a number of bulk formulas +% evaluated by Fung et al (1984), Rev. Geophys. Space Physics, +% 22,177-193. +% +% INPUT: Ts - sea surface temperature [C] +% Ta - air temperature [C] +% rh - relative humidity [%] +% Fc - cloudiness correction factor F(C) (=1 for clear sky) +% (see cloudcor.m for more information) +% bulk_f - bulk formulas to be used +% (see Fung et al, 1984) for details). +% Options are: +% 'brunt' +% 'berliand' - probably the best one (default) +% 'clark' +% 'hastenrath' +% 'efimova' +% 'bunker' +% 'anderson' +% 'swinbank' - does not use rh +% +% OUTPUT: lwest - net downward longwave flux [W/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/31/98: version 1.1 (contributed by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% load constants and defaults +as_consts; + +if nargin<5, + bulk_f=bulkf_default; +end; + +% compute vapour pressure from relative humidity (using formulas +% from Gill, 1982) +ew=satvap(Ta); +rw=eps_air*ew./(P_default-ew); +r=(rh/100).*rw; +e_a=r*P_default./(eps_air+r); + +% convert to K +ta=Ta+CtoK; +ts=Ts+CtoK; + +% signs for flux INTO ocean +switch bulk_f, + case {'brunt',1} + lwest = -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc; + + case {'berliand',2} + lwest = -emiss_lw*sigmaSB.*ta.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc ... + - 4*emiss_lw*sigmaSB.*ta.^3.*(Ts - Ta); + + case {'clark',3} + lwest = -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc ... + - 4*emiss_lw*sigmaSB.*ts.^3.*(Ts - Ta); + + case {'hastenrath',4} + q_a = e_a*eps_air./P_default; + lwest = -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.056*sqrt(q_a)).*Fc ... + - 4*emiss_lw*sigmaSB.*ts.^3.*(Ts - Ta); + + case {'efimova',5} + lwest = -emiss_lw*sigmaSB.*ta.^4.*(0.254 - 0.00495*(e_a)).*Fc ; + + case {'bunker',6} + lwest = -0.022*emiss_lw*sigmaSB.*ta.^4.*(11.7 - 0.23*(e_a)).*Fc ... + - 4*emiss_lw*sigmaSB.*ta.^3.*(Ts - Ta); + + case {'anderson',5} + lwest = -emiss_lw*sigmaSB.*(ts.^4 - ta.^4.*(0.74+0.0049*(e_a))).*Fc ; + + case {'swinbank',5} + lwest = -emiss_lw*sigmaSB.*(ts.^4 - 9.36e-6*ta.^6).*Fc ; + + otherwise + error(['Unrecognized bulk formula specified: ' bulk_f]); +end; + diff --git a/cdnlp.m b/cdnlp.m new file mode 100644 index 0000000..27d3b30 --- /dev/null +++ b/cdnlp.m @@ -0,0 +1,34 @@ +function [cd,u10]=cdnlp(sp,z) +% CDNLP: computes neutral drag coefficient following Large&Pond (1981). +% [cd,u10]=CDNLP(sp,z) computes the neutral drag coefficient and wind +% speed at 10m given the wind speed at height z following Large and +% Pond (1981), J. Phys. Oceanog., 11, 324-336. +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% +% OUTPUT: cd - neutral drag coefficient at 10m +% u10 - wind speed at 10m [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +as_consts; % define physical constants +a=log(z./10)/kappa; % log-layer correction factor +tol=.001; % tolerance for iteration [m/s] + +u10o=zeros(size(sp)); +cd=1.15e-3*ones(size(sp)); +u10=sp./(1+a.*sqrt(cd)); + +ii=abs(u10-u10o)>tol; +while any(ii(:)), + u10o=u10; + cd=(4.9e-4+6.5e-5*u10o); % compute cd(u10) + cd(u10o<10.15385)=1.15e-3; + u10=sp./(1+a.*sqrt(cd)); % next iteration + ii=abs(u10-u10o)>tol; % keep going until iteration converges +end; diff --git a/cdntc.m b/cdntc.m new file mode 100644 index 0000000..0403d5e --- /dev/null +++ b/cdntc.m @@ -0,0 +1,52 @@ +function [cd,u10]=cdntc(sp,z,Ta) +% CTDTC: computes the neutral drag coefficient following Smith (1988). +% [cd,u10]=CDNTC(sp,z,Ta) computes the neutral drag coefficient and +% wind speed at 10m given the wind speed and air temperature at height z +% following Smith (1988), J. Geophys. Res., 93, 311-326. +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% Ta - air temperature (optional) [C] +% +% OUTPUT: cd - neutral drag coefficient at 10m +% u10 - wind speed at 10m [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +if nargin==2, + Ta=Ta_default; +end; + +% iteration endpoint +tol=.00001; + +visc=viscair(Ta); + +% remove any sp==0 to prevent division by zero +i=find(sp==0); +sp(i)=.1.*ones(length(i),1); + +% initial guess +ustaro=zeros(size(sp)); +ustarn=.036.*sp; + +% iterate to find z0 and ustar +ii=abs(ustarn-ustaro)>tol; +while any(ii(:)), + ustaro=ustarn; + z0=Charnock_alpha.*ustaro.^2./g + R_roughness*visc./ustaro; + ustarn=sp.*(kappa./log(z./z0)); + ii=abs(ustarn-ustaro)>tol; +end + +sqrcd=kappa./log((10)./z0); +cd=sqrcd.^2; + +u10=ustarn./sqrcd; diff --git a/cdnve.m b/cdnve.m new file mode 100644 index 0000000..2c0b42a --- /dev/null +++ b/cdnve.m @@ -0,0 +1,43 @@ +function [cd,u10]=cdnve(sp,z) +% CDNVE: computes neutral drag coefficient following Vera (1983). +% [cd,u10]=CDNVE(sp,z) computes the neutral drag coefficient and wind +% speed at 10m given the wind speed at height z. Uses the expression +% for friction velocity derived by E. Vera (1983) and published as +% eqn. 8 in Large, Morzel, and Crawford (1995), J. Phys. Oceanog., 25, +% 2959-2971. Range of fit to data is 1 to 25 m/s. +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% +% OUTPUT: cd - neutral drag coefficient at 10m +% u10 - wind speed at 10m [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (modified by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% constants in fit for drag coefficient +A=2.717e-3; +B=0.142e-3; +C=0.0764e-3; + +as_consts; % other constants +a=log(z./10)/kappa; % log-layer correction factor +tol=.001; % tolerance for iteration (m/s) + +u10o=zeros(size(sp))+.1; % don't start iteration at 0 to prevent blowups. +cd=(A./u10o + B + C*u10o); + +u10=sp./(1+a.*sqrt(cd)); + +ii=abs(u10-u10o)>tol; +while any(ii(:)), + u10o=u10; + cd=(A./u10o + B + C*u10o); + u10=sp./(1+a.*sqrt(cd)); % next iteration + ii=abs(u10-u10o)>tol; % keep going until iteration converges +end; + + diff --git a/cloudcor.m b/cloudcor.m new file mode 100644 index 0000000..896b9b0 --- /dev/null +++ b/cloudcor.m @@ -0,0 +1,70 @@ +function Fc=cloudcor(C,optns,lat) +% CLOUDCOR: computes cloud correction factor for bulk long-wave flux. +% Fc=CLOUDCOR(C,optns,lat) computes the cloud correction factor F(C) +% as a function of the cloud fraction C for bulk long-wave flux formulae. +% In general, these are functions of the form +% 1 - a_n*C^n +% Since the coefficients and powers depend a lot on the dominant cloud +% type which may vary from region to region and season to season, it is +% not clear which parametrization is best (see Fung et al (1984), +% Rev. of Geophys. and Space Phys., 22, 177-193). +% +% The particular parametrization used here depends on the second input +% variable, for which no default is given to emphasize the fact that you +% really need to understand what you are doing here! +% +% optns = [a1 a2] = use a correction factor of [1-a1*C-a2*C^2]. +% +% There are several "built-in" formulae (from Fung et al) that all have +% a latitude-dependence of some kind. +% +% optns = 'clarke',lat = Clarke (1974) corrections for abs(lat)<50. +% = 'bunker',lat = Bunker (1976) corrections for N Atlantic. +% +% INPUT: C - cloud fraction +% optns - see above for details +% lat - latitude [deg] - required for "built-in" formulae only +% +% OUTPUT: Fc - correction factor used as input to BLWHF + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/12/98: version 1.1 (contributed by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if isstr(optns), + + switch optns(1), + case 'c', + a1=0; + if abs(lat)>55, a2=NaN; + elseif abs(lat)>45, a2=0.73; + elseif abs(lat)>35, a2=0.69; + elseif abs(lat)>25, a2=0.64; + elseif abs(lat)>15, a2=0.60; + elseif abs(lat)> 7, a2=0.56; + elseif abs(lat)> 2, a2=0.53; + else a2=0.51; + end; + + case 'b', + a2=0; + if abs(lat)>75, a1=0.84; + elseif abs(lat)>65, a1=0.80; + elseif abs(lat)>55, a1=0.76; + elseif abs(lat)>45, a1=0.72; + elseif abs(lat)>35, a1=0.68; + elseif abs(lat)>25, a1=0.63; + elseif abs(lat)>15, a1=0.59; + elseif abs(lat)>7, a1=0.52; + else a1=0.50; + end; + + otherwise + error('Unrecognized option'); + end; +else + a1=optns(1);a2=optns(2); +end; + +Fc = 1 - a1*C - a2.*C.^2; diff --git a/clskswr.m b/clskswr.m new file mode 100644 index 0000000..4e64f7c --- /dev/null +++ b/clskswr.m @@ -0,0 +1,66 @@ +function sradav = clskswr(yd,lat) +% CLSKSWR: computes clear sky insolation following Seckel&Beaudry (1973). +% sradav = CLSKSWR(yd,lat) computes average clear sky solar insolation +% based on the Seckel and Beaudry (1973) formula presented in Reed (1977), +% J. Phys. Ocean., 7, 482-485. Assumes the year is not a leap year. +% +% INPUT: yd - yearday (e.g., Jan 10th is 10) +% lat - latitude [deg] +% +% OUTPUT: sradav - clear sky mean daily insolation [W/m^2] + +% NOTE: The output appears to be very similar to what you would get by +% averaging soradna.m output over a day and then assuming an +% atmospheric transmission of 0.7 (with differences of order 10% +% for latitudes below 40N, and increasing to 30% in winter at 60N). +% In absolute terms the agreement is to within +/-25 W/m^2. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/28/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +radlat = pi*lat/180; + +if length(lat)==1, + lat=lat+zeros(size(yd)); +end; + +% check if yd is negative +ind=find(yd<0); +yd(ind)=yd(ind)+365; +% truncate to integer yearday for use with formula +yd=fix(yd); + +phi = (yd-21)*360/365; +phi = pi*phi/180; + +sradav=zeros(size(yd))+zeros(size(lat))+NaN; + +ii= lat>=-20 & lat<40; +if any(ii(:)), + a0 = -15.82 + 326.87*cos(radlat(ii)); + a1 = 9.63 + 192.44*cos(radlat(ii)+pi/2); + b1 = -3.27 + 108.70*sin(radlat(ii)); + a2 = -0.64 + 7.80*sin(2*(radlat(ii)-pi/4)); + b2 = -0.50 + 14.42*cos(2*(radlat(ii)-5*pi/180)); + sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii)); +end + +ii= lat>=40 & lat<=60;; +if any(ii(:)), + l2=lat(ii).^2; + a0 = 342.61 - 1.97*lat(ii) - 0.018*l2; + a1 = 52.08 - 5.86*lat(ii) + 0.043*l2; + b1 = -4.80 + 2.46*lat(ii) -0.017*l2; + a2 = 1.08 - 0.47*lat(ii) + 0.011*l2; + b2 = -38.79 + 2.43*lat(ii) - 0.034*l2; + sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii)); +end + +if any(lat>60 | lat<-20) + warning('Formula only works for latitudes 20S-60N, see help text for further help') +end + + diff --git a/contents.m b/contents.m new file mode 100644 index 0000000..7286b22 --- /dev/null +++ b/contents.m @@ -0,0 +1,56 @@ +% Contents of AIR_SEA TOOLBOX, Version 2.0 +% 8/9/99 +% +% +% air_dens.m 966 08-09-99 12:14p air_dens.m +% albedo.m 1,009 08-09-99 11:49a albedo.m +% albedot1.mat 7,757 07-21-98 5:59p albedot1.mat +% as_consts.m 3,699 08-09-99 11:54a as_consts.m +% binave.m 1,209 08-09-99 11:50a binave.m +% blwhf.m 2,775 08-05-99 4:54p blwhf.m +% cdnlp.m 1,197 08-05-99 5:13p cdnlp.m +% cdntc.m 1,356 08-05-99 5:18p cdntc.m +% cdnve.m 1,418 08-05-99 5:26p cdnve.m +% cloudcor.m 2,418 08-05-99 5:24p cloudcor.m +% clskswr.m 2,244 08-05-99 5:25p clskswr.m +% cool_skin.m 7,331 08-09-99 11:46a cool_skin.m +% delq.m 1,059 08-09-99 11:34a delq.m +% ep.m 852 08-09-99 12:12p ep.m +% greg2.m 1,511 08-09-99 11:35a greg2.m +% hfbulktc.m 12,232 08-09-99 11:39a hfbulktc.m +% hms2h.m 714 08-09-99 11:39a hms2h.m +% julianmd.m 1,534 08-09-99 11:47a julianmd.m +% lfhf.m 1,626 08-09-99 11:47a lwhf.m +% no_skin_out.txt 13,953 08-05-99 3:41p no_skin_out.txt +% omegalmc.m 1,021 08-09-99 11:40a omegalmc.m +% qsat.m 1,329 08-09-99 11:52a qsat.m +% rain_flux.m 2,569 08-05-99 3:47p rain_flux.m +% readme.m 8,114 08-09-99 2:06p readme.m +% reedcf.m 3,807 08-05-99 3:50p reedcf.m +% relhumid.m 1,519 08-09-99 11:52a relhumid.m +% rhadj.m 707 08-05-99 3:53p rhadj.m +% satvap.m 785 08-09-99 11:42a satvap.m +% s2hms.m 478 08-05-99 3:53p s2hms.m +% slhftc.m 2,031 08-09-99 11:56a slhftc.m +% soradna1.m 4,486 08-09-99 11:42a soradna1.m +% spshftlp.m 993 08-09-99 11:43a spshftlp.m +% spshfttc.m 1,118 08-09-99 11:43a spshfttc.m +% spshftve.m 991 08-09-99 11:43a spshftve.m +% stresstc.m 1,061 08-09-99 11:57a stresstc.m +% stresslp.m 990 08-09-99 11:56a stresslp.m +% stressve.m 974 08-09-99 11:57a stressve.m +% sunrise.m 1,821 08-09-99 11:58a sunrise.m +% swhf.m 1,418 08-05-99 4:19p swhf.m +% t_hfbulktc.m 1,704 08-05-99 4:45p t_hfbulktc.m +% test2_5b.dat 11,699 08-05-99 4:20p test2_5b.dat +% vapor.m 1,019 08-09-99 11:44a vapor.m +% viscair.m 569 08-05-99 4:33p viscair.m +% wavdist2.m 876 08-05-99 4:37p wavdist2.m +% wavdist1.m 1,432 08-05-99 4:35p wavdist1.m +% wavdist.m 1,517 08-05-99 4:39p wavedist.m +% wdnotes.m 1,745 08-05-99 4:42p wdnotes.m +% yearday.m 1,162 08-05-99 4:43p yearday.m +% yes_skin_out.txt 13,958 08-05-99 4:43p yes_skin_out.txt + + + diff --git a/cool_skin.m b/cool_skin.m new file mode 100644 index 0000000..9981ceb --- /dev/null +++ b/cool_skin.m @@ -0,0 +1,187 @@ +function [delta,Dter,Dqer] = cool_skin(sal,Tw,rhoa,cpa,Pa, ... + U_star,T_star,Q_star, ... + dlw,dsw,nsw,delta,g,Rgas, ... + CtoK,Qsat_coeff) +% COOL_SKIN: compute the cool-skin parameters. +% COOL_SKIN computes the cool-skin parameters. This code follows +% the fortran program bulk_v25b.f. For more details, see the cool-skin +% and warm layer paper by Fairall et al (1996), JGR, 101, 1295-1308. +% All input variables should be vectors (either row or column), except +% Rgas, CtoK, Qsat_coeff, and g, which can be scalars. Uses some +% functions from CSIRO SEAWATER TOOLBOX. +% +% INPUT: sal - salinity [psu (PSS-78)] +% Tw - water surface temperature [C] +% rhoa - air density [kg/m^3] +% cpa - specific heat capacity of air [J/kg/C] +% Pa - air pressure [mb] +% U_star - friction velocity including gustiness [m/s] +% T_star - temperature scale [C] +% Q_star - humidity scale [kg/kg] +% dlw - downwelling (INTO water) longwave radiation [W/m^2] +% dsw - measured insolation [W/m^2] +% nsw - net shortwave radiation INTO water [W/m^2] +% delta - cool-skin layer thickness [m] +% g - gravitational constant [m/s^2] +% Rgas - gas constant for dry air [J/kg/K] +% CtoK - conversion factor for deg C to K +% Qsat_coeff - saturation specific humidity coefficient +% +% OUTPUT: delta - cool-skin layer thickness [m] +% Dter - cool-skin temperature difference [C]; positive if +% surface is cooler than bulk (presently no warm skin +% permitted by model) +% Dqer - cool-skin specific humidity difference [kg/kg] +% +% USAGE: [delta,Dter,Dqer] = cool_skin(sal,Tw,rhoa,cpa,Pa, ... +% U_star,T_star,Q_star, ... +% dlw,dsw,nsw,delta,g,Rgas, ... +% CtoK,Qsat_coeff) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 4/9/99: version 1.2 (contributed by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% -> column vectors +sal = sal(:); Tw = Tw(:); rhoa = rhoa(:); cpa = cpa(:); Pa = Pa(:); +U_star = U_star(:); T_star = T_star(:); Q_star = Q_star(:); +dlw = dlw(:); nsw = nsw(:); delta = delta(:); Rgas = Rgas(:); +CtoK = CtoK(:); Qsat_coeff = Qsat_coeff(:); g = g(:); + + +size_data = size(Tw); + +alpha = sw_alpha(sal,Tw,0); % thermal expansion coeff [1/C] +beta_sal = sw_beta(sal,Tw,0); % saline contraction coeff [1/psu] +cpw = sw_cp(sal,Tw,0); % specific heat capacity [J/kg/C] +rhow = sw_dens0(sal,Tw); % density at atmospheric press [kg/m^3] +viscw = sw_visc(sal,Tw,0); % kinematic viscosity of sea-water [m^2/s] +tcondw = sw_tcond(sal,Tw,0); % thermal conductivity of sea-water [W/m/K] + +% the following are values used for COARE +% alpha = 2.1e-5*(Tw+3.2).^0.79;% as used for COARE data +% beta_sal = 0.026./sal; % as used for COARE data +% cpw = 4000*ones(size(Tw)); % as used for COARE data +% rhow = 1022*ones(size(Tw)); % as used for COARE data +% viscw = 1e-6*ones(size(Tw)); % as used for COARE data +% tcondw = 0.6*ones(size(Tw)); % as used for COARE data + +% latent heat of water +Le = (2.501-0.00237*Tw)*10^6; + +% saturation specific humidity; +Qs = Qsat_coeff.*qsat(Tw,Pa); + +% a big constant +bigc = (16.*g.*cpw.*(rhow.*viscw).^3)./(tcondw.^2.*rhoa.^2); + +% constant for correction of dq; slope of sat. vap. +wetc = 0.622.*Le.*Qs./(Rgas.*(Tw+CtoK).^2); + +% compute fluxes out of the ocean (i.e., up = positive) +hsb = - rhoa.*cpa.*U_star.*T_star; +hlb = - rhoa.*Le.*U_star.*Q_star; + +% net longwave (positive up) +nlw = - lwhf(Tw,dlw,dsw); + +% total heat flux out of the water surface +qout = nlw + hsb + hlb; + +% compute deltaSc = fc*Sns, see sec. 2.4 (p. 1297-1298) in cool-skin paper +% 3 choices; comment out those that are not used! +deltaSc = zeros(size_data); +ipos_nsw = find(nsw > 0); +deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),1).*nsw(ipos_nsw); % Paulson and Simpson (1981) +% deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),2).*nsw(ipos_nsw); % COARE approx. to Paulson +% deltaSc(ipos_nsw) = f_c(delta(ipos_nsw),3).*nsw(ipos_nsw); % Hasse (1971) + +qcol = qout - deltaSc; + +% initialize +alphaQb = zeros(size_data); +lamda = zeros(size_data); +Dter = zeros(size_data); + +ipos_qcol = find(qcol > 0); + +% eqn. 17 in cool-skin paper +alphaQb(ipos_qcol) = alpha(ipos_qcol).*qcol(ipos_qcol) + ... + sal(ipos_qcol).*beta_sal(ipos_qcol) ... + .*hlb(ipos_qcol).*cpw(ipos_qcol)./Le(ipos_qcol); + +% eqn. 14 in cool-skin paper +lamda(ipos_qcol) = 6./(1+(bigc(ipos_qcol).*alphaQb(ipos_qcol) ... + ./U_star(ipos_qcol).^4).^0.75).^(1/3); + +% eqn. 12 in cool-skin paper +delta(ipos_qcol) = lamda(ipos_qcol).*viscw(ipos_qcol) ... + ./(sqrt(rhoa(ipos_qcol)./rhow(ipos_qcol)) ... + .*U_star(ipos_qcol)); + +% eqn. 13 in cool-skin paper +Dter(ipos_qcol) = qcol(ipos_qcol).*delta(ipos_qcol)./tcondw(ipos_qcol); + +Dqer = wetc.*Dter; + + + + +function fc = f_c(delta,option) +% F_C: computes the absorption coefficient fc. +% fc=F_C(delta,option) computes the absorption coefficient fc. +% +% INPUT: delta - thickness of cool-skin layer [m] +% option - 1 use Paulson and Simpson (1981) data for seawater; +% See also p. 1298 of Fairall et al (1996) JGR, 101, +% cool-skin and warm-layer paper. +% +% 2 use approximation to Paulson as given in +% p. 1298 of Fairall et al (1996) JGR, 101, cool-skin +% and warm-layer paper. +% +% 3 use fc = const. = 0.19, as suggested by Hasse (1971). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/5/99: version 1.2 (contributed by AA) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if option == 1 + % Use Paulson and Simpson data + + % Wavelength bands for the coefficients [um] + % 1) 0.2-0.6 + % 2) 0.6-0.9 + % 3) 0.9-1.2 + % 4) 1.2-1.5 + % 5) 1.5-1.8 + % 6) 1.8-2.1 + % 7) 2.1-2.4 + % 8) 2.4-2.7 + % 9) 2.7-3.0 + + % F_i is the amplitude + F_i = [0.237 0.360 0.179 0.087 0.080 0.0246 0.025 0.007 0.0004]; + F_i1 = repmat(F_i,length(delta),1); + + % Gam_i is the absorption length [m] + Gam_i = [34.8 2.27 3.15e-2 5.48e-3 8.32e-4 1.26e-4 3.13e-4 7.82e-5 1.44e-5]; + Gam_i1 = repmat(Gam_i,length(delta),1); + + delta1 = repmat(delta,1,length(Gam_i)); + + % fc is the absorption in the cool-skin layer of thickness delta + fc = sum(F_i1.*(1-(Gam_i1./delta1).*(1-exp(-delta1./Gam_i1))), 2); + +elseif option == 2 + % use COARE approximation to Paulson and Simpson data + + fc = 0.137+11.*delta-(6.6e-5./delta).*(1-exp(-delta/8e-4)); + +elseif option == 3 + % use Hasse simple approximation + + fc = 0.19; + +end diff --git a/delq.m b/delq.m new file mode 100644 index 0000000..31c08f2 --- /dev/null +++ b/delq.m @@ -0,0 +1,26 @@ +function dq=delq(Ts,Ta,rh) +% DELQ: computes air-sea specific humidity difference. +% dq=DELQ(Ts,Ta,rh) computes the specific humidity (kg/kg) difference +% between the air (as determined by relative humidty rh and air +% temperature Ta measurements) and the sea surface (where q is +% assumed to be at 98% satuation at the sea surface temperature Ts). +% DELQ uses QSAT based on Tetens' formula for saturation vapor +% pressure from Buck (1981), J. App. Meteor., 1527-1532. The +% dependence of QSAT on pressure is small (<0.5%) and has been +% removed using a mean pressure of 1020 mb. +% +% INPUT: Ts - sea surface temperature [C] +% Ta - air temperature [C] +% rh - relative humidity [%] +% +% OUTPUT: dq - air-sea specific humidity difference [kg/kg] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 4/10/98: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +dq=0.01.*rh.*qsat(Ta) - 0.98.*qsat(Ts); + + diff --git a/ep.m b/ep.m new file mode 100644 index 0000000..53a9267 --- /dev/null +++ b/ep.m @@ -0,0 +1,29 @@ +function [E,P]=ep(rfd,Qlat) +% EP: computes evap and precip accumulation +% [E,P]=EP(rfd,Qlat) computes precipitation and evaporation +% accumulation from rainfall rate rfd and latent heat +% flux Qlat. Assumes hourly input. +% +% INPUT: rfd - precip rate [mm/min] +% Qlat - latent heat flux [W/m^2] +% +% OUTPUT: P - precip accumulation [m] +% E - evaporation accumulation [m] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% compute P +P=cumsum(rfd)./1000; % convert mm to m +dt=60; % for hourly data, dt = 60 min +P=P.*dt; + +% compute E +Le=2.5e6; % heat of vaporization (W/m^2) +pw=1025; % density of seawater (kg/m^3) at 32 psu, 10 degC, 0 db +dt=3600; % seconds per hour +E=cumsum(Qlat).*dt./(Le.*pw); + + + diff --git a/greg2.m b/greg2.m new file mode 100644 index 0000000..4282f5f --- /dev/null +++ b/greg2.m @@ -0,0 +1,53 @@ + function [gtime]=greg2(yd,yr) +% GREG2: converts decimal yearday to standard Gregorian time. +% [gtime]=GREG2(yd,yr) converts decimal yearday to corresponding +% Gregorian calendar dates. In this convention, Julian day 2440000 +% begins at 0000 UT May 23 1968. +% +% INPUT: yd - decimal yearday (e.g., 0000 UT Jan 1 is 0.0) +% yr - year (e.g., 1995) +% +% OUTPUT: gtime is a six component Gregorian time vector +% gtime=[year mo da hr mi sec] +% +% Example: [1995 01 01 12 00 00] = greg2(0.5, 1995) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/10/97: version 1.0 +% 4/7/99: version 1.2 (simplified by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +js = julianmd(yr,01,01,00); +julian = js + yd; +julian=julian+5.e-9; % kludge to prevent roundoff error on seconds + +% if you want Julian Days to start at noon... +% h=rem(julian,1)*24+12; +% i=(h >= 24); +% julian(i)=julian(i)+1; +% h(i)=h(i)-24; Otherwise,.... + +secs=rem(julian,1)*24*3600; + +j = floor(julian) - 1721119; +in = 4*j -1; +y = floor(in/146097); +j = in - 146097*y; +in = floor(j/4); +in = 4*in +3; +j = floor(in/1461); +d = floor(((in - 1461*j) +4)/4); +in = 5*d -3; +m = floor(in/153); +d = floor(((in - 153*m) +5)/5); +y = y*100 +j; +mo=m-9; +yr=y+1; +i=(m<10); +mo(i)=m(i)+3; +yr(i)=y(i); +[hour,min,sec]=s2hms(secs); +gtime=[yr(:) mo(:) d(:) hour(:) min(:) sec(:)]; + + diff --git a/hfbulktc.m b/hfbulktc.m new file mode 100644 index 0000000..7598fc4 --- /dev/null +++ b/hfbulktc.m @@ -0,0 +1,360 @@ +function A=hfbulktc(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw) +% 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]; + +% Code follows Edson and Fairall TOGA COARE code (version 2.0), modified +% to include Rogers' weighting factor for unstable conditions. Code does +% include gustiness, and assumes that the marine boundary layer height is +% known and constant over time for simiplicity. zr/L is limited to +% be <=3.0 to ensure that the code converges to nonzero stress and heat +% flux values for strongly stable conditions. The bulk Richardson number +% is computed between the sea surface and zr as a diagnostic about whether +% turbulent boundary layer theory is applicable. Code does not include +% warm layer effects to modify Ts. See Fairall et al (1996), J. Geophys. +% Res., 101, 3747-3764, for description of full TOGA COARE code and +% comparison with data. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/19/98: version 1.1 (rewritten by RP to remove inconsistencies in +% virtual and real temperatures, improve loop structure, +% correct gustiness component of stress computation) +% 4/9/99: version 1.2 (rewritten by AA to clarify some variable names +% and include cool-skin effect and Webb correction to latent +% heat flux added to output matrix) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +M=length(ur); + +% change to column vectors +ur = ur(:); +Ta = Ta(:); +rh = rh(:); +Ts = Ts(:); +Pa = Pa(:); +zr = zr(:); +zt = zt(:); +zq = zq(:); + +% create vectors for rh, Pa, zr, zt, and zq, if scalars are input +if length(rh)==1 & M>1 + rh=rh*ones(M,1); +end +if length(Pa)==1 & M>1 + Pa=Pa*ones(M,1); +end +if length(zr)==1 & M>1 + zr=zr*ones(M,1); +end +if length(zt)==1 & M>1 + zt=zt*ones(M,1); +end +if length(zq)==1 & M>1 + zq=zq*ones(M,1); +end + +if nargin > 8 + % optional cool-skin stuff + sal = sal(:); + dlw = dlw(:); + dsw = dsw(:); + nsw = nsw(:); + % create vector for sal if scalar is input + if length(sal)==1 & M>1 + sal=sal*ones(M,1); + end +end + +% initialize various constants +as_consts; + +tol=.001; % tolerance on Re changes to make sure soln has converged. + +onethird=1./3; +o61=1/eps_air-1; % 0.61 (moisture correction for temperature) + +visc=viscair(Ta); % viscosity +Qsats=Qsat_coeff*qsat(Ts,Pa); % saturation specific humidity; the Qsat_coeff + % value is set in routine as_consts.m +Q=(0.01.*rh).*qsat(Ta,Pa); % specific humidity of air [kg/kg] +T =Ta+CtoK; % convert to K +Tv=T.*(1 + o61*Q); % air virtual temperature +rho=(100*Pa)./(gas_const_R*Tv); % air density +Dt=(Ta+0.0098.*zt)-Ts; % adiabatic temperature difference +Dq=Q-Qsats; % humidity difference + +% compute initial neutral scaling coefficients +S=sqrt(ur.^2 + min_gustiness.^2); +cdnhf=sqrt(cdntc(S,zr,Ta)); % Smith's neutral cd as first guess + +z0t=7.5*10^(-5); +ctnhf=kappa./log(zt./z0t); + +z0q=z0t; +cqnhf=kappa./log(zq./z0q); + +U_star = cdnhf.*S; % (includes gustiness) +T_star = ctnhf.*Dt; % +Q_star = cqnhf.*Dq; % + +Dter = 0; +Dqer = 0; +if nargin > 8 +% initial cool-skin thickness guess + delta = 0.001*ones(size(Ts)); +end + +Reu=0;Ret=0;Req=0; +% begin iteration loop to compute best U_star, T_star, and Q_star +for iter1=1:80; + + ReuO=Reu; RetO=Ret; ReqO=Req; % Save old values + + % Compute Monin-Obukov length (NB - definition given as eqn (7) + % of Fairall et al (1996) probably wrong, following, e.g. + % Godfrey and Bellars (1991), JGR, 96, 22043-22048 and original code) + bs=g*(T_star.*(1 + o61*Q) + o61*T.*Q_star)./Tv; + L=(U_star.^2)./(kappa*bs); + % set upper limit on zr/L = 3.0 to force convergence under + % very stable conditions. Assume that zr, zt and zq comparable. + index_limit = (L0); + L(index_limit)=zr(index_limit)/3; + + zetu=zr./L; % nondimensionalized heights + zett=zt./L; + zetq=zq./L; + + % surface roughness + z0=(Charnock_alpha/g).*U_star.^2 + R_roughness.*visc./U_star; + + % compute U_star correction for non-neutral conditions + cdnhf=kappa./(log(zr./z0)-psiutc(zetu)); + U_star=cdnhf.*S; + + Reu=z0.*U_star./visc; % roughness Reynolds # + [Ret,Req]=LKB(Reu); % compute other roughness Reynolds #s + + % compute t and q roughness scales from roughness R#s + z0t=visc.*Ret./U_star; + z0q=visc.*Req./U_star; + + % compute new transfer coefficients at measurement heights + cthf=kappa./(log(zt./z0t)-psittc(zett)); + cqhf=kappa./(log(zq./z0q)-psittc(zetq)); + + % compute new values of T_star, Q_star + T_star=cthf.*(Dt + Dter); + Q_star=cqhf.*(Dq + Dqer); + + % estimate new gustiness + Ws=U_star.*(-CVB_depth./(kappa*L)).^onethird; + wg=min_gustiness*ones(M,1); + j=find(zetu<0); % convection in unstable conditions only + wg(j)=max(min_gustiness,beta_conv.*Ws(j)); % set minimum gustiness + S=sqrt(ur.^2 + wg.^2); + + if nargin > 8 + % compute cool-skin parameters + [delta,Dter,Dqer] = cool_skin(sal,Ts-Dter,rho,cp,Pa, ... + U_star,T_star,Q_star, ... + dlw,dsw,nsw,delta,g,gas_const_R, ... + CtoK,Qsat_coeff); + end + +end % end of iteration loop + +ii= abs(Reu-ReuO)>tol | abs(Ret-RetO)>tol | abs(Req-ReqO)>tol; +if any(ii), + disp(['Algorithm did not converge for ' int2str(sum(ii)) ' values. Indices are:']); + disp(find(ii)'); + warning('Not converged!'); +end; + + +% compute latent heat +Le=(2.501-0.00237*(Ts-Dter))*10^6; + +% compute fluxes into ocean +Hs=rho.*cp.*U_star.*T_star; +Hl=rho.*Le.*U_star.*Q_star; + +% compute transfer coefficients at measurement heights +CD=(U_star./S).^2; +CT=U_star.*T_star./(S.*(Dt + Dter)); % Stanton number +CQ=U_star.*Q_star./(S.*(Dq + Dqer)); % Dalton number + +% to compute mean stress, we don't want to include the effects +% of gustiness which average out (in a vector sense). +stress=rho.*CD.*S.*ur; + +% compute bulk Richardson number (as a diagnostic) - the "T" +% is probably not quite right - assumes T \ approx. Ts (good enough though) +RI=g.*zr.*((Dt + Dter) + o61*T.*(Dq + Dqer))./(Tv.*S.^2); + +% compute Webb correction to latent heat flux into ocean +W = 1.61.*U_star.*Q_star + (1 + 1.61.*Q).*U_star.*T_star./T; % eqn. 21 +Hl_webb = rho.*Le.*W.*Q; % eqn. 22, Fairall et al. (1996), JGR, 101, p3751. + +% output array +if nargin > 8 + % output additional cool-skin parameters + A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter]; +else + % otherwise + A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI]; +end + + + +function y=psiutc(zet) +% PSIUTC: computes velocity profile function following TOGA/COARE. +% y=PSIUTC(zet) computes the turbulent velocity profile function given +% zet = (z/L), L the Monin-Obukoff length scale, following Edson and +% Fairall TOGA COARE code (version 2.0) as modified to include Rogers' +% weighting factor to combine the Dyer and free convection forms for +% unstable conditions. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 +% 8/5/99: version 1.2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +c13=1.0./3.0; +sq3=sqrt(3.0); + +% stable conditions +y=-4.7*zet; + +% unstable conditions +j=find(zet<0); +zneg=zet(j); + +% nearly stable (standard functions) + x=(1-16.0.*zneg).^0.25; + y1=2.0.*log((1+x)./2) + log((1+x.^2)./2) -2.*atan(x) + pi/2; + +% free convective limit + x=(1-12.87*zneg).^c13; + y2=1.5*log((x.^2+x+1)./3) - sq3*atan((2.*x+1)/sq3) + pi/sq3; + +% weighted sum of the two + F=1.0./(1+zneg.^2); + y(j)=F.*y1+(1-F).*y2; + + + + +function y=psittc(zet) +% PSITTC: computes potential temperature profile following TOGA/COARE. +% y=PSITTC(zet) computes the turbulent potential temperature profile +% function given zet = (z/L), L the Monin-Obukoff length scale, following +% Edson and Fairall TOGA COARE code (version 2.0), as modified to use +% Rogers' weighting factor to combine the Dyer and free convective +% forms for unstable conditions. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 +% 8/5/99: version 1.2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +c13=1.0./3.0; +sq3=sqrt(3.0); + +% stable conditions +y=-4.7.*zet; + +% unstable conditions +j=find(zet<0); +zneg=zet(j); + +% nearly stable (standard functions) + x=(1-16.0.*zneg).^0.25; + y1=2.0*log((1+x.^2)./2); + +% free convective limit + x=(1-12.87*zneg).^c13; + y2=1.5.*log((x.^2+x+1)./3.0) - sq3.*atan((2.*x+1)./sq3) + pi./sq3; + +% weighted sum of the two + F=1.0./(1+zneg.^2); + y(j)=F.*y1 + (1-F).*y2; + + + + + +function [Ret,Req]=LKB(Reu); +% LKB: computes rougness Reynolds numbers for temperature and humidity +% [Ret,Req]=LKB(Reu) computes the roughness Reynolds for temperature +% and humidity following Liu, Katsaros and Businger (1979), J. Atmos. +% Sci., 36, 1722-1735. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 +% 8/5/99: version 1.2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + Ret=.177*ones(size(Reu)); + Req=.292*ones(size(Reu)); +j=find(Reu>.11 & Reu<=.825); + Ret(j)=1.376.*Reu(j).^0.929; + Req(j)=1.808.*Reu(j).^0.826; +j=find(Reu>.825 & Reu<=3); + Ret(j)=1.026./Reu(j).^0.599; + Req(j)=1.393./Reu(j).^0.528; +j=find(Reu>3 & Reu<=10); + Ret(j)=1.625./Reu(j).^1.018; + Req(j)=1.956./Reu(j).^0.870; +j=find(Reu>10 & Reu<=30); + Ret(j)=4.661./Reu(j).^1.475; + Req(j)=4.994./Reu(j).^1.297; +j=find(Reu>30); + Ret(j)=34.904./Reu(j).^2.067; + Req(j)=30.790./Reu(j).^1.845; diff --git a/hms2h.m b/hms2h.m new file mode 100644 index 0000000..3a988c2 --- /dev/null +++ b/hms2h.m @@ -0,0 +1,29 @@ +function [hours]=hms2h(h,m,s); +% HMS2H: converts hours, minutes, and seconds to hours. +% [hours]=HMS2H(h,m,s) converts hours, minutes, and seconds to hours +% +% INPUT: h - integer hours +% m - minutes +% s - secs +% +% OUTPUT: hours - decimal hours +% +% Usage: [hours]=hms2h(h,m,s); or [hours]=hms2h(hhmmss); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin== 1, + hms=h; + h=floor(hms/10000); + ms=hms-h*10000; + m=floor(ms/100); + s=ms-m*100; + hours=h+m/60+s/3600; +else + hours=h+(m+s/60)/60; +end + + diff --git a/julianmd.m b/julianmd.m new file mode 100644 index 0000000..cebdfe4 --- /dev/null +++ b/julianmd.m @@ -0,0 +1,47 @@ +function [j]=julianmd(y,m,d,h) +% JULIANMD: converts Gregorian calendar time to decimal Julian day. +% [j]=JULIANMD(y,m,d,h) converts Gregorian calendar dates to corresponding +% Julian day numbers. Although the formal definition holds that Julian +% days start and end at noon, here Julian days start and end at midnight. +% In this convention, Julian day 2440000 began at 0000 UT, May 23, 1968. +% +% INPUT: d - day (1-31) component of Gregorian date +% m - month (1-12) component +% y - year (e.g., 1979) component +% h - decimal hours (assumed 0 if absent) +% +% OUTPUT: j - decimal Julian day number (e.g., 0000 UT Jan 1 is 0.0) +% +% Usage: [j]=julianmd(y,m,d,h) (inputs scalars or matrices) +% or +% [j]=julianmd([y m d hr min sec]) (nx6 matrix input) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin==3, + h=0.; +elseif nargin==1, + h=hms2h(y(:,4),y(:,5),y(:,6)); + d=y(:,3); + m=y(:,2); + y=y(:,1); +end +mo=m+9; +yr=y-1; +i=(m>2); +mo(i)=m(i)-3; +yr(i)=y(i); +c = floor(yr/100); +yr = yr - c*100; +j = floor((146097*c)/4) + floor((1461*yr)/4) + ... + floor((153*mo +2)/5) +d +1721119; + +% if you want Julian days to start and end at noon, +% replace the following line with: +% j=j+(h-12)/24; + +j=j+h/24; + diff --git a/lwhf.m b/lwhf.m new file mode 100644 index 0000000..225bd03 --- /dev/null +++ b/lwhf.m @@ -0,0 +1,47 @@ +function qlw=lwhf(Ts,dlw,dsw) +% LWHF: computes net longwave heat flux following Dickey et al (1994). +% qlw=LWHF(Ts,dlw) computes the net longwave heat flux into the ocean. +% Following Dickey et al (1994), J. Atmos. Oceanic Tech., 11, 1057-1078, +% the incident longwave flux can be corrected for sensor heating due to +% insolation if you are using Epply or Kipp & Zonen CG1 pyrgeometers. +% In this case, use qlw=LWHF(Ts,dlw,dsw). Epply is the default +% pyrgeometer; change code for the Kipp & Zonen instrument. +% +% INPUT: Ts - sea surface temperature [C] +% dlw - (measured) downward longwave flux [W/m^2] +% dsw - (measured) insolation [W/m^2] (needed for Eppley +% or Kipp & Zonen pyrgeometers) +% +% OUTPUT: qlw - net longwave heat flux [W/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/19/98: version 1.1 (revised for non-Epply pyrgeometers by RP) +% 4/9/99: version 1.2 (included Kipp & Zonen CG1 pyrgeometers by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +% convert degC to degK +ts=Ts+CtoK; + +% correct dlw for sensor heating by insolation +if nargin==3, + % this line is for Epply pyrgeometers + dlwc=dlw-0.036.*dsw; + + % this line is for Kipp & Zonen CG1 pyrgeometers + % (the offset is specified as 25 W/m^2 at 1000 W/m^2) + % dlwc=dlw-0.025.*dsw; +else + dlwc=dlw; +end; + +% compute upward gray-body longwave flux +lwup=-emiss_lw.*sigmaSB.*(ts.^4); + +% compute net flux +qlw=lwup + emiss_lw.*dlwc; + diff --git a/no_skin_out.txt b/no_skin_out.txt new file mode 100644 index 0000000..babb67c --- /dev/null +++ b/no_skin_out.txt @@ -0,0 +1,119 @@ +% Output from Fortran program bulk_v25b.f, with both cool-skin and warm-layer +% turned OFF, and using input file test2_5b.dat with COARE data +% index xtime hsb hlb tub ts HF EF TAU T0 Webb RainF rain dt_cool dt_warm tk_pwp + 1 921125132100. 6. 111. .030 29.15 9. 125. .030 29.15 5. 0. .0 .00 .00 19.00 + 2 921125141200. 5. 96. .023 29.15 8. 108. .023 29.15 4. 0. .0 .00 .00 19.00 + 3 921125150300. 5. 98. .024 29.15 7. 110. .025 29.15 4. 0. .0 .00 .00 19.00 + 4 921125155500. 5. 110. .029 29.15 8. 124. .029 29.15 5. 0. .0 .00 .00 19.00 + 5 921125164600. 4. 92. .019 29.16 7. 104. .019 29.16 4. 0. .0 .00 .00 19.00 + 6 921125173800. 9. 123. .036 29.16 13. 138. .035 29.16 6. 0. .0 .00 .00 19.00 + 7 921125182900. 4. 61. .009 29.16 6. 70. .009 29.16 3. 0. .0 .00 .00 19.00 + 8 921125192000. 4. 89. .019 29.15 6. 99. .019 29.15 4. 0. .0 .00 .00 19.00 + 9 921125201200. 4. 113. .030 29.14 7. 125. .030 29.14 5. 0. .0 .00 .00 19.00 + 10 921125210500. 5. 84. .021 29.13 7. 93. .020 29.13 4. 0. .0 .00 .00 19.00 + 11 921125221900. 5. 113. .028 29.14 7. 121. .028 29.14 4. 0. .0 .00 .00 19.00 + 12 921125232700. 6. 124. .036 29.14 7. 128. .036 29.14 5. 0. .0 .00 .00 19.00 + 13 921126001900. 6. 106. .026 29.20 6. 108. .025 29.20 4. 0. .0 .00 .00 19.00 + 14 921126011000. 5. 90. .022 29.27 6. 93. .021 29.27 4. 0. .0 .00 .00 19.00 + 15 921126022200. 4. 123. .030 29.26 4. 121. .030 29.26 4. 0. .0 .00 .00 19.00 + 16 921126031300. 5. 119. .030 29.42 6. 123. .030 29.42 4. 0. .0 .00 .00 19.00 + 17 921126041600. 4. 104. .025 29.59 6. 110. .024 29.59 4. 0. .0 .00 .00 19.00 + 18 921126050700. 3. 94. .025 29.59 5. 103. .025 29.59 4. 0. .0 .00 .00 19.00 + 19 921126055900. 2. 90. .024 29.50 5. 103. .025 29.50 4. 0. .0 .00 .00 19.00 + 20 921126065500. 2. 96. .022 29.37 4. 107. .022 29.37 4. 0. .0 .00 .00 19.00 + 21 921126074700. 1. 93. .020 29.37 4. 106. .020 29.37 4. 0. .0 .00 .00 19.00 + 22 921126083800. 1. 78. .014 29.35 3. 86. .014 29.35 3. 0. .0 .00 .00 19.00 + 23 921126093000. 1. 56. .007 29.39 2. 63. .007 29.39 2. 0. .0 .00 .00 19.00 + 24 921126102100. 1. 50. .006 29.45 3. 58. .006 29.45 2. 0. .0 .00 .00 19.00 + 25 921126111200. 1. 51. .006 29.42 3. 58. .006 29.42 2. 0. .0 .00 .00 19.00 + 26 921126120800. 1. 60. .009 29.42 3. 70. .009 29.42 3. 0. .0 .00 .00 19.00 + 27 921126130000. 2. 71. .012 29.35 4. 79. .012 29.35 3. 0. .0 .00 .00 19.00 + 28 921126145300. 3. 58. .011 29.41 5. 68. .010 29.41 3. 0. .0 .00 .00 19.00 + 29 921126154500. 4. 65. .012 29.31 6. 74. .013 29.31 3. 0. .0 .00 .00 19.00 + 30 921126163600. 4. 68. .013 29.29 6. 76. .013 29.29 3. 0. .0 .00 .00 19.00 + 31 921126172700. 2. 57. .009 29.27 4. 66. .009 29.27 3. 0. .0 .00 .00 19.00 + 32 921126181900. 2. 47. .006 29.26 4. 55. .006 29.26 2. 0. .0 .00 .00 19.00 + 33 921126191000. 2. 36. .003 29.27 3. 41. .002 29.27 2. 0. .0 .00 .00 19.00 + 34 921126200100. 1. 30. .002 29.27 2. 35. .002 29.27 1. 0. .0 .00 .00 19.00 + 35 921126213400. 2. 31. .002 29.29 2. 34. .001 29.29 1. 0. .0 .00 .00 19.00 + 36 921126222600. 5. 55. .007 29.38 7. 61. .007 29.38 3. 0. .0 .00 .00 19.00 + 37 921126231700. 39. 182. .093 29.31 43. 198. .093 29.31 12. 26. 4.8 .00 .00 19.00 + 38 921127000900. 30. 175. .064 29.33 33. 186. .063 29.33 10. 0. .0 .00 .00 19.00 + 39 921127010000. 19. 160. .048 29.33 22. 171. .048 29.33 8. 0. .0 .00 .00 19.00 + 40 921127015100. 14. 129. .032 29.33 16. 138. .032 29.33 6. 0. .0 .00 .00 19.00 + 41 921127024300. 12. 118. .029 29.36 15. 130. .028 29.36 6. 0. .0 .00 .00 19.00 + 42 921127033400. 15. 129. .037 29.25 18. 140. .038 29.25 7. 0. .0 .00 .00 19.00 + 43 921127042500. 35. 161. .070 29.23 39. 177. .070 29.23 11. 50. 9.4 .00 .00 19.00 + 44 921127051700. 32. 154. .048 29.23 36. 169. .048 29.23 10. 9. 1.6 .00 .00 19.00 + 45 921127060800. 54. 227. .154 29.24 60. 248. .155 29.24 16. 9. 1.5 .00 .00 19.00 + 46 921127070000. 29. 137. .040 29.22 33. 152. .040 29.22 9. 0. .0 .00 .00 19.00 + 47 921127075100. 16. 121. .027 29.22 20. 135. .027 29.22 6. 0. .0 .00 .00 19.00 + 48 921127084200. 10. 112. .020 29.16 13. 123. .020 29.16 5. 0. .0 .00 .00 19.00 + 49 921127093400. 7. 102. .016 29.15 9. 112. .016 29.15 4. 0. .0 .00 .00 19.00 + 50 921127102500. 6. 113. .021 29.14 9. 126. .021 29.14 4. 0. .0 .00 .00 19.00 + 51 921127113200. 7. 132. .034 29.14 11. 147. .034 29.14 6. 0. .0 .00 .00 19.00 + 52 921127122600. 8. 147. .048 29.12 12. 163. .047 29.12 6. 0. .0 .00 .00 19.00 + 53 921127134800. 7. 105. .033 29.11 10. 118. .032 29.11 5. 0. .0 .00 .00 19.00 + 54 921127143900. 7. 110. .036 29.10 10. 123. .036 29.10 5. 0. .0 .00 .00 19.00 + 55 921127160400. 7. 126. .039 29.13 10. 140. .039 29.13 5. 0. .0 .00 .00 19.00 + 56 921127170000. 5. 102. .026 29.14 8. 115. .026 29.14 4. 0. .0 .00 .00 19.00 + 57 921127175200. 4. 89. .020 29.12 7. 100. .020 29.12 4. 0. .0 .00 .00 19.00 + 58 921127184300. 4. 75. .013 29.12 6. 85. .013 29.12 3. 0. .0 .00 .00 19.00 + 59 921127193500. 3. 69. .012 29.12 5. 78. .012 29.12 3. 0. .0 .00 .00 19.00 + 60 921127202600. 2. 54. .007 29.12 4. 61. .007 29.12 2. 0. .0 .00 .00 19.00 + 61 921127211700. 3. 54. .006 29.12 4. 57. .006 29.12 2. 0. .0 .00 .00 19.00 + 62 921127220900. 3. 52. .005 29.13 4. 53. .004 29.13 2. 0. .0 .00 .00 19.00 + 63 921127230000. 5. 53. .004 29.14 4. 49. .004 29.14 2. 0. .0 .00 .00 19.00 + 64 921127235200. 7. 56. .003 29.16 4. 44. .003 29.16 2. 0. .0 .00 .00 19.00 + 65 921128004300. 8. 63. .003 29.14 4. 44. .003 29.14 2. 0. .0 .00 .00 19.00 + 66 921128022100. 9. 74. .005 29.17 4. 50. .004 29.17 2. 0. .0 .00 .00 19.00 + 67 921128042600. 6. 66. .004 29.27 2. 44. .003 29.27 2. 0. .0 .00 .00 19.00 + 68 921128043600. 5. 69. .004 29.47 2. 50. .003 29.47 2. 0. .0 .00 .00 19.00 + 69 921128052800. 3. 51. .002 29.23 0. 36. .001 29.23 1. 0. .0 .00 .00 19.00 + 70 921128061900. 2. 44. .002 29.44 1. 34. .001 29.44 1. 0. .0 .00 .00 19.00 + 71 921128075200. 2. 50. .003 29.27 1. 42. .003 29.27 1. 0. .0 .00 .00 19.00 + 72 921128084300. 2. 45. .003 29.30 1. 39. .002 29.30 1. 0. .0 .00 .00 19.00 + 73 921128093500. 4. 71. .010 29.30 3. 68. .009 29.30 2. 0. .0 .00 .00 19.00 + 74 921128102600. 2. 56. .005 29.35 2. 59. .005 29.35 2. 0. .0 .00 .00 19.00 + 75 921128111800. 1. 42. .002 29.24 1. 42. .002 29.24 1. 0. .0 .00 .00 19.00 + 76 921128120900. 3. 61. .011 29.30 5. 67. .011 29.30 3. 0. .0 .00 .00 19.00 + 77 921128130100. 3. 56. .008 29.56 5. 63. .007 29.56 3. 0. .0 .00 .00 19.00 + 78 921128135200. 2. 47. .004 29.56 4. 55. .004 29.56 2. 0. .0 .00 .00 19.00 + 79 921128144400. 2. 44. .004 29.29 3. 48. .004 29.29 2. 0. .0 .00 .00 19.00 + 80 921128153500. 2. 52. .006 29.26 3. 58. .006 29.26 2. 0. .0 .00 .00 19.00 + 81 921128162600. 2. 57. .007 29.26 4. 66. .007 29.26 2. 0. .0 .00 .00 19.00 + 82 921128171800. 2. 56. .008 29.19 4. 64. .008 29.19 2. 0. .0 .00 .00 19.00 + 83 921128180900. 2. 48. .006 29.13 3. 55. .006 29.13 2. 0. .0 .00 .00 19.00 + 84 921128190100. 2. 58. .008 29.18 4. 65. .008 29.18 2. 0. .0 .00 .00 19.00 + 85 921128195200. 2. 55. .007 29.17 3. 62. .007 29.17 2. 0. .0 .00 .00 19.00 + 86 921128204400. 2. 49. .005 29.17 3. 54. .004 29.17 2. 0. .0 .00 .00 19.00 + 87 921128213500. 3. 59. .006 29.40 4. 62. .006 29.40 2. 0. .0 .00 .00 19.00 + 88 921128222600. 6. 46. .003 29.52 6. 46. .003 29.52 2. 0. .0 .00 .00 19.00 + 89 921128231800. 6. 48. .003 29.52 5. 44. .003 29.52 2. 0. .0 .00 .00 19.00 + 90 921129000900. 7. 47. .002 29.58 5. 36. .001 29.58 2. 0. .0 .00 .00 19.00 + 91 921129010100. 10. 90. .011 29.49 6. 73. .010 29.49 3. 0. .0 .00 .00 19.00 + 92 921129015200. 9. 84. .011 29.41 6. 74. .009 29.41 3. 0. .0 .00 .00 19.00 + 93 921129024400. 6. 69. .006 29.41 3. 55. .005 29.41 2. 0. .0 .00 .00 19.00 + 94 921129033500. 5. 77. .008 29.37 2. 64. .007 29.37 2. 0. .0 .00 .00 19.00 + 95 921129042600. 4. 83. .009 29.23 1. 67. .008 29.23 2. 0. .0 .00 .00 19.00 + 96 921129051800. 3. 70. .007 29.23 1. 59. .006 29.23 2. 0. .0 .00 .00 19.00 + 97 921129060900. 4. 54. .005 29.32 3. 50. .004 29.32 2. 0. .0 .00 .00 19.00 + 98 921129070100. 17. 112. .021 29.26 14. 102. .020 29.26 5. 32. 6.5 .00 .00 19.00 + 99 921129075200. 13. 75. .009 29.46 14. 80. .008 29.46 4. 41. 6.6 .00 .00 19.00 + 100 921129091700. 7. 76. .008 29.34 7. 76. .008 29.34 3. 0. .0 .00 .00 19.00 + 101 921129100800. 6. 93. .018 29.33 6. 93. .017 29.33 4. 0. .0 .00 .00 19.00 + 102 921129110000. 3. 74. .014 29.23 4. 77. .014 29.23 3. 0. .0 .00 .00 19.00 + 103 921129115500. 2. 60. .009 29.23 3. 65. .009 29.23 2. 0. .0 .00 .00 19.00 + 104 921129124700. 2. 60. .010 29.18 4. 66. .010 29.18 2. 0. .0 .00 .00 19.00 + 105 921129133800. 3. 54. .007 29.40 4. 61. .007 29.40 2. 0. .0 .00 .00 19.00 + 106 921129142900. 3. 57. .007 29.40 4. 61. .007 29.40 2. 0. .0 .00 .00 19.00 + 107 921129153100. 4. 66. .009 29.63 6. 72. .009 29.63 3. 0. .0 .00 .00 19.00 + 108 921129163800. 5. 79. .013 29.57 6. 86. .013 29.57 3. 0. .0 .00 .00 19.00 + 109 921129173000. 5. 84. .016 29.49 7. 91. .016 29.49 4. 0. .0 .00 .00 19.00 + 110 921129182100. 4. 77. .014 29.42 6. 85. .013 29.42 3. 0. .0 .00 .00 19.00 + 111 921129191300. 3. 71. .012 29.42 6. 81. .012 29.42 3. 0. .0 .00 .00 19.00 + 112 921129200400. 5. 78. .014 29.35 7. 85. .014 29.35 3. 0. .0 .00 .00 19.00 + 113 921129205500. 3. 73. .012 29.35 5. 80. .012 29.35 3. 0. .0 .00 .00 19.00 + 114 921129214700. 4. 67. .009 29.35 5. 70. .009 29.35 3. 0. .0 .00 .00 19.00 + 115 921129223800. 5. 71. .009 29.37 5. 70. .009 29.37 3. 0. .0 .00 .00 19.00 + 116 921129233000. 7. 77. .009 29.31 5. 70. .009 29.31 3. 0. .0 .00 .00 19.00 diff --git a/omegalmc.m b/omegalmc.m new file mode 100644 index 0000000..0efe6cf --- /dev/null +++ b/omegalmc.m @@ -0,0 +1,30 @@ +function y=omegalmc(x) +% OMEGALMC: estimates wind log profile correction due to surface waves. +% y=OMEGALMC(x) computes the log profile correction function due to wind +% distortion associated with surface waves. Input is x=za/Hw, where za +% is the measurement height and Hw is the dominant surface wave height. +% Functional form is simplified (analytic) version of empirical omega +% curves shown in Fig. 9b of Large, Morzel, and Crawford (1995), J. Phys. +% Oceanog., 25, 2959-2971, with the wave-induced roughness length xr=0.15. +% Assumes x is a vector with all elements greater than zr. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +xr=0.15; +ylimit=-6; +y=ylimit.*ones(size(x)); +i=find(x<3.2967); +% polynomial fit +a=-2.6; +p1=-0.0199; +p2=0.0144; +p3=0.7660; +p4=0.0654; +x2=x(i).^2;x3=x2.*x(i); +y(i)=a.*log(x(i)./xr)+p1.*x3+p2.*x2+p3.*x(i)+p4; + + + diff --git a/qsat.m b/qsat.m new file mode 100644 index 0000000..d42cb39 --- /dev/null +++ b/qsat.m @@ -0,0 +1,35 @@ + function q=qsat(Ta,Pa) +% QSAT: computes specific humidity at saturation. +% q=QSAT(Ta) computes the specific humidity (kg/kg) at satuation at +% air temperature Ta (deg C). Dependence on air pressure, Pa, is small, +% but is included as an optional input. +% +% INPUT: Ta - air temperature [C] +% Pa - (optional) pressure [mb] +% +% OUTPUT: q - saturation specific humidity [kg/kg] + +% Version 1.0 used Tetens' formula for saturation vapor pressure +% from Buck (1981), J. App. Meteor., 1527-1532. This version +% follows the saturation specific humidity computation in the COARE +% Fortran code v2.5b. This results in an increase of ~5% in +% latent heat flux compared to the calculation with version 1.0. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 4/7/99: version 1.2 (revised as above by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if nargin==1, + as_consts; + Pa=P_default; % pressure in mb +end; + +% original code +% a=(1.004.*6.112*0.6220)./Pa; +% q=a.*exp((17.502.*Ta)./(240.97+Ta)) + +% as in Fortran code v2.5b for COARE +ew = 6.1121*(1.0007+3.46e-6*Pa).*exp((17.502*Ta)./(240.97+Ta)); % in mb +q = 0.62197*(ew./(Pa-0.378*ew)); % mb -> kg/kg diff --git a/rain_flux.m b/rain_flux.m new file mode 100644 index 0000000..c5e820c --- /dev/null +++ b/rain_flux.m @@ -0,0 +1,59 @@ +function [tau_rain,heat_rain] = rain_flux(Ta,Pa,rh,rain,Ts,sal,u,zu) +% RAIN_FLUX: computes heat and momentum flux due to rain. +% RAIN_FLUX computes heat flux and momentum flux due to rain. This +% code follows the Fortran program bulk_v25b.f. For more details, +% see Fairall et al. (1996), JGR, 101, 3751-3752. +% +% INPUT: Ta - air temperature [C] +% Pa - air pressure [mb] +% rh - relative humidity [%] +% rain - rain rate [mm/hr] +% Ts - sea surface temperature [C] +% sal - salinity [psu (PSS-78)] +% u - wind speed [m/s] +% zu - wind measurement height [m] +% +% OUTPUT: tau_rain - momentum flux of rainfall [N/m^2] +% heat_rain - heat flux of rainfall (OUT of ocean) [W/m^2] +% +% USAGE: [tau_rain,heat_rain] = RAIN_FLUX(Ta,Pa,rh,rain,Ts,sal,u,zu) + +% NOTE: All input variables should be vectors (either row or column), zu +% may also be a fixed scalar. Output variables are column vectors. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 4/3/99: version 1.2 (contributed by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% -> column vectors +Ta = Ta(:); Pa = Pa(:); rh = rh(:); rain = rain(:); Ts = Ts(:); +sal = sal(:); u = u(:); zu = zu(:); + +% get constants +as_consts; +cpa = cp; + + +o61 = 1/eps_air - 1; % ~0.61 (moisture correction for temp.) +Qa = 0.01*rh.*qsat(Ta,Pa); % specific humidity of air [kg/kg] +T = Ta + CtoK; % C -> K +Tv = T.*(1 + o61*Qa); % air virtual temperature +rhoa = (100*Pa)./(gas_const_R*Tv); % air density +Le = (2.501-0.00237*Ts)*1e6; % latent heat of vaporization at Ts +Qs = Qsat_coeff*qsat(Ts,Pa); % saturation specific humidity + + +% compute heat flux of rainfall OUT of ocean +dwat = 2.11e-5*(T./CtoK).^1.94; % water vapour diffusivity +dtmp = (1+3.309e-3*Ta-1.44e-6*Ta.^2)*0.02411./(rhoa.*cpa); % heat diffusivity +dqs_dt = Qa.*Le./(gas_const_R*T.^2); % Clausius-Clapeyron +alfac = 1./(1+0.622*(dqs_dt.*Le.*dwat)./(cpa.*dtmp)); % wet bulb factor +cpw = sw_cp(sal,Ts,0); % heat capacity of sea water + +heat_rain = rain.*alfac.*cpw.*((Ts-Ta)+(Qs-Qa).*Le./cpa)/3600; + +% compute momentum flux of rainfall +[cd10,u10] = cdntc(u,zu,Ta);% use Smith's formula to compute wind speed at 10m +tau_rain = rain.*u10/3600; + diff --git a/readme.m b/readme.m new file mode 100644 index 0000000..66e93af --- /dev/null +++ b/readme.m @@ -0,0 +1,180 @@ +% AIR_SEA: Introduction to the AIR_SEA TOOLBOX + +% AIR_SEA TOOLBOX (version 2.0: 8/9/99) +% +% 1) Introduction: Welcome to the AIR_SEA toolbox, a collection of MATLAB +% programs (m-files) which can be used to compute surface wind stress and +% heat flux components from buoy and shipboard atmospheric and near-surface +% oceanographic time series measurements. All m-files include a header +% which describes the mfile's function, input and output variables, and +% key references when important. They have been written for use with +% MATLAB 5. +% +% 2) Conventions: While not required for many m-files, it is generally +% assumed that the input time series of measured variables are hourly +% averaged column or row vectors and the other input variables are scalars, +% all expressed in MKS units. Two time conventions are used: a) decimal +% Julian yearday where 0000 UT Jan 1 equals 0.0, and b) calender yearday +% where Jan 1 equals 1. The choice of which convention is used is +% internally consistent between m-files. +% +% 3) Programs used to compute heat flux components: +% +% Shortwave flux: +% +% SWHF: computes net sw flux into ocean and albedo. Uses +% SORADNA1 and ALBEDO to compute solar altitude, no-sky +% insolation and albedo. +% +% SORADNA1: computes no-sky insolation and solar altitude at +% a given time and location. +% +% ALBEDO: computes ocean albedo following Payne (1972). +% +% Longwave flux: +% +% LWHF: computes net lw flux into ocean when downward +% lw radiation is measured, using Dickey et al (1994). +% +% BLWHF: computes net lw flux into the ocean when downward +% lw radiation is NOT measured. Uses SATVAP. Requires +% as input a cloudiness correction factor from CLOUDCOR. +% +% CLOUDCOR: cloudiness correction factor used in bulk formulae, +% based on estimated Cloud Fraction, which is either observed +% directly or estimated, using, e.g., REEDCF. +% +% REEDCF: computes daily average Cloud Fraction using formula of +% Reed (1977), who relates daily average cloudiness to the observed +% reduction in solar insolation from clear-sky values. +% +% Sensible and latent fluxes: +% +% HFBULKTC: uses a simplified version of Fairall et al (1996) +% TOGA/COARE code to compute sensible and latent heat +% fluxes into ocean. Present version includes a) Rogers' +% weighting factor for unstable conditions, b) the effects +% of gustiness, c) a constant marine boundary layer height, +% d) a limit of zr/L <=3.0 to ensure that the code converges +% to nonzero stress and heat flux values for strongly stable +% conditions, e) cool-skin effect, and f) Webb correction for +% latent heat flux. NOTE: both cool-skin and Webb correction +% are optional, and user must decide if they want these used, +% e.g., in SLHFTC. Warm layer effects are not included in this +% version. Uses VISCAIR and QSAT to compute air viscosity +% and saturation specific humidity, CDNTC the neutral drag +% coefficient, and PSIUTC and PSITTC to adjust the different +% transfer coefficients to the measurement heights for a +% given stability. Also returns related variables. +% +% SLHFTC: includes ocean surface current and HFBULKTC to comput +% sensible and latent heat fluxes into ocean. +% +% RAIN_FLUX: computes heat flux and momentum flux due to rain. +% +% 4) Programs relating wind speed, height, and surface stress. +% +% Neutral conditions: +% +% CDNLP: computes neutral Cd, 10m wind following Large and Pond (1981). +% CDNTC: computes neutral Cd, 10m wind following Smith (1988). +% CDNVE: computes neutral Cd, 10m wind following Vera (1983). +% +% STRESSLP: computes the neutral wind stress using Large and Pond. +% STRESSTC: computes the neutral wind stress following Smith. +% STRESSVE: computes the neutral wind stress using Vera. +% +% SPSHFTLP: computes winds at another height using Large&Pond drag. +% SPSHFTTC: computes winds at another height using Smith drag. +% SPSHFTVE: computes winds at another height using Vera drag. +% +% Non-neutral conditions: +% +% HFBULKTC: uses simplified version of Fairall et al (1996) +% TOGA/COARE code to compute surface wind stress amplitude, +% (Uses Monin-Obukov similarity theory with surface rougness using +% Charnock approach, like Smith (1988)). +% +% SLHFTC: includes ocean surface current and HFBULKTC to +% compute surface wind stress vector as well as scalar parameters. +% +% 5) Programs used to estimate wave effects on the measured wind speed: +% +% WAVEDIST: estimate true wind speed at 10-m height. +% WAVEDIS1: estimate true wind speed at measurement height. +% WAVEDIS2: plots wave effects at measurement height vs. wave height. +% OMEGALMC: estimates wave effect on wind log profile. +% CDNVE: computes neutral drag coefficient following Vera (1983). +% +% See WDNOTES for additional information. +% +% 6) Other useful programs: +% +% AS_CONSTS: contains various constants used in the toolbox. +% +% DAVEALB: computes daily mean albedo. +% SUNRISE: computes GMT time of sunrise and sunset (uses SORADNA1). +% +% GREG2: converts decimal yearday into Julian calendar day. +% JULIANMD: converts Gregorian calendar dates to decimal Julian day +% for days beginning at midnight UT +% YEARDAY: converts calender month and day into yearday. +% +% DELQ: air-sea specific humidity difference. +% EP: net precipitation and evaporation accumulation. +% QSAT: saturation specific humidity. +% RELHUMID: relative humidity from wet/dry bulb thermometers. +% RHADJ: adjusts RH for values above 100. +% SATVAP: saturation vapour pressure. +% VAPOR: heat of evaporation. +% VISCAIR: viscosity of air at a given temperature. +% COOL_SKIN: computes cool-skin parameters. +% T_HFBULKTC: tests HFBULKTC with COARE data. +% +% 7) See CONTENTS for listing of all m-files in this toolbox. +% +% 8) History: +% +% Version 1.0: +% +% The initial assembly of this toolbox was a collaborative effort +% by Bob Beardsley (WHOI), Ed Dever (SIO), Steve Lentz (WHOI), Jim +% Edson (WHOI), and Dick Payne (WHOI), with additional input from +% Steve Anderson (WHOI), Jay Austin (WHOI), Chris Fairall (NOAA), +% Carl Friehe (UCI), Bill Large (NCAR), Dave Rogers (SIO), Rich +% Signell (USGS), and Bob Weller (WHOI). Their input was very useful. +% +% Version 1.1: +% +% Rich Pawlowicz (UBC) then converted the original version 1.0 +% (written for MATLAB 4) into a much improved version 1.1 (optimized +% for MATLAB 5) which included major coding improvements, the addition +% of some new m-files, and some corrections of existing m-files. +% +% Version 1.2: +% +% Ayal Anis (U. Dalhousie) then modified HFBULKTC to include the +% Fairall et al (1996) cool-skin effect and Webb correction to the +% latent heat flux, plus added files to test the code with COARE +% data. Ayal and R. Onken (NATO) also contributed several other files. +% +% Version 2.0: +% +% Bob Beardsley has added several m-files and made simple changes +% to the various m-files to standardize the format and documentation. +% +% 9) Comments, Suggestions, and Improvements +% +% Please contact Bob Beardsley at rbeardsley@whoi.edu with questions +% and comments, especially concerning bugs (and their possible fixes), +% ideas for additional m-files, plus any m-files which you want to +% contribute to this toolbox. Your help in improving this toolbox will +% be greatly appreciated. +% +% As new or/or improved m-files are developed for this toolbox, they +% will be added to the AIR_SEA toolbox folder located at the SEA-MAT +% Web site (crusty.usgs.gov/sea-mat/). SEA-MAT is a collection of +% MATLAB mfiles for the display and analysis of oceanographic data. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + diff --git a/reedcf.m b/reedcf.m new file mode 100644 index 0000000..2616ceb --- /dev/null +++ b/reedcf.m @@ -0,0 +1,136 @@ +function c = reedcf(yd,lat,dsw); +% REEDCF: computes daily mean cloud cover following Reed (1977). +% c = REEDCF(yd,lat,dsw) computes daily averaged cloud cover c from +% yearday, latitude, and measured insolation following Reed (1977), J. Phys. +% Oceanog., 7, 482-485. Assumes hourly input series are either both +% column or both row vectors of equal length. c is output as a boxcar +% function with c constant over a 24 hr period. +% +% INPUT: yd - yearday (e.g., Jan 10 is yd=10) +% lat - latitude [deg] +% dsw - (measured) insolation [W/m^2] +% +% OUTPUT: c - daily averaged cloud cover + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% determine if input is either row or column vectors +[m,n] = size(dsw); + +% first daily average swr and yd +davesw = binave(dsw,24); +daveyd = binave(yd,24); + +% estimate daily averaged cloud factor +c=cloudfac(daveyd,lat,davesw); + +% cloudfac always outputs a column vector, if initial input is +% row vector, convert cf to row vector +if m==1 + c = c'; +end + +% convert daily averaged cloudfactor to boxcar function +if n== 1 + c = c*ones(1,24); + c = reshape(c',prod(size(c)),1); + +elseif m==1 + c = ones(24,1)*c; + c = reshape(c,1,prod(size(c))); +end +c(length(dsw)+1:end)=[]; + + +function cf=cloudfac(yd,lat,davesw) +% CLOUDFAC: computes daily mean cloud factor following Reed (1977). +% cf=CLOUDFAC(yd,lat,davesw) computes the daily average cloud factor +% based on Reed (1977), J. Phys. Ocean., 7, 482-485. Assumes year is +% not a leap year. If yd is negative, then yd=yd+365. +% +% INPUT: yd - yearday (e.g., Jan 10th is 10) +% lat - latitude [deg] +% davesw - average measured insolation [W/m^2] +% +% OUTPUT: cf - daily average cloud factor + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% convert to column vectors +[N,M]=size(davesw); +if N==1 + davesw=davesw'; + yd=yd'; +end + +% check if yd is negative +ind=find(yd<0); +yd(ind)=yd(ind)+365; +% truncate to integer yearday for use with formula +yd=fix(yd); + +% determine noon solar altitude +nsa=nsunang(yd,lat); + +% compute ratio of observed to clear sky insolation +cssw=clskswr(yd,lat); +QdQ=davesw./cssw; + +% correct for measurement error + +ind1=find(QdQ>1); +QdQ(ind1)=ones(length(ind1),1); + +% compute cloud factor +cf=(1-QdQ+.0019.*nsa)./.62; + +% truncate limits + +ind2=find(cf>1); +cf(ind2)=ones(length(ind2),1); +ind3=find(cf<0.3); +cf(ind3)=zeros(length(ind3),1); + +% reconvert if necessary +if N==1 + cf=cf'; +end + + +function nsa=nsunang(yd,lat) +% NSUNANG: computes noon solar altitude angle. +% nsa=NSUNANG(yd,lat) computes the noonday solar altitude angle as a +% function of yearday and latitude using Smithsonian table (see +% Reed (1976), NOAA Tech Memo., ERL PMEL-8, 20 pp., or Reed (1977), +% J. Phys. Ocean., 7, 482-485). +% +% INPUT: yd - yearday (e.g., Jan 10 is 10) +% lat - latitude [deg] +% +% OUTPUT: nsa - noonday solar altitude [deg] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ensure that yd is positive +ind=find(yd<0); +yd(ind)=yd(ind)+365; +yd=fix(yd); + +% compute noon solar altitude +k=pi./180; +t=2.*pi.*(yd./365); +d=.397+3.630*sin(t)-22.98.*cos(t)+.040.*sin(2.*t)-.388.*cos(2.*t)... + +.075.*sin(3.*t)-.160.*cos(3.*t); +dl=k.*d; +ll=k.*lat; +z=sin(ll).*sin(dl)+cos(ll).*cos(dl); +nsa=asin(z)./k; diff --git a/relhumid.m b/relhumid.m new file mode 100644 index 0000000..3c0cf17 --- /dev/null +++ b/relhumid.m @@ -0,0 +1,54 @@ +function rh=relhumid(Td,Tw,P,p_typ) +% RELHUMID: finds relative humidity from wet/dry thermometer readings. +% rh=relhumid(Td,Tw,Pa,type) computes the relative humidity from +% wt and dry-bulb temperature measurements using the psychrometric eqn. +% and constants from Sargent (1980), Meteorol. Mag. 109, 238-246. The +% latter two inputs are optional. +% +% INPUTS : Td - dry bulb thermometer [C] +% Tw - wet thermometer [C] +% Pa - air pressure (optional) [mb] +% type - 'assman' for Assman-type forced ventilation +% 'screen' for standard screen (natural ventilation) +% +% OUTPUT: rh - relative humidity [%] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 (contributed by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +as_consts; + +if nargin==2, + P=P_default; + p_typ=psych_default; +elseif nargin==3, + if isstr(P), + p_typ=P; + P=P_default; + else + p_typ=psych_default; + end; +end; + +% psychrometric coefficient + +switch p_typ, + case 'screen', + A=0.000799; % natural screens + case 'assman', + A = 0.000667; % Assmann-type with forced ventilation + otherwise + error(['unknown psychrometer type: ' p_typ]); +end; + +% compute saturation vapour pressure for both temps. +ed=satvap(Td,P); +ewp=satvap(Tw,P); + +% The psychrometric eqn! +e = ewp - A*P.*(Td-Tw); % ambient vapour pressure + +rh= e./ed * 100; + diff --git a/rhadj.m b/rhadj.m new file mode 100644 index 0000000..616b4bd --- /dev/null +++ b/rhadj.m @@ -0,0 +1,22 @@ +function rha = rhadj(rh,rhmax) +% RHADJ: rescales RH to have a maximum of 100%. +% RHADJ(rh,rhmax) rescales RH so that the maximum values do not +% exceed 100%. Assumes values between 93% and rhmax should be +% rescaled to 93 - 100% (the calibration curves of RH sensors usually +% become nonlinear above ~ 90%, and may peak above 100% in this range +% above ~ 90% where their calibration becomes unstable.) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 4/10/98: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +rhl=93; + +rhn=rh; +a=(100-rhl)./(rhmax-rhl); +drh=rh-rhl; +j=find(drh>0); +rhn(j)=rhl+a.*drh(j); + +rha=rhn; \ No newline at end of file diff --git a/s2hms.m b/s2hms.m new file mode 100644 index 0000000..d6ae96f --- /dev/null +++ b/s2hms.m @@ -0,0 +1,15 @@ +function [hr,min,sec]=s2hms(secs) +% S2HMS: converts seconds to interger hour, minute, and seconds. +% [hr,min,sec]=S2HMS(secs) converts seconds to integer hour, minute, +% and seconds. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/11/96: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +sec=round(secs); +hr=floor(sec./3600); +min=floor(rem(sec,3600)./60); +sec=round(rem(sec,60)); + diff --git a/satvap.m b/satvap.m new file mode 100644 index 0000000..9ec1a27 --- /dev/null +++ b/satvap.m @@ -0,0 +1,30 @@ +function ew=satvap(Ta,P) +% SATVAP: computes saturation vapor pressure +% q=satvap(Ta) computes the vapor pressure at satuation at air +% temperature Ta (deg C). From Gill (1982), Atmos-Ocean Dynamics, 606. +% +% INPUT: Ta- air temperature [C] +% p - pressure (optional) [mb] +% +% OUTPUT: q - saturation vapour pressure [mb] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/27/98: version 1.1 (corrected by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +as_consts; + +if nargin<2, + P=P_default; +end; +if nargin==0, + Ta=Ta_default; +end; + +ew=power(10,((0.7859+0.03477*Ta)./(1+0.00412*Ta))); + +fw=1 + 1e-6*P.*(4.5+0.0006*Ta.^2); + +ew=fw.*ew; diff --git a/slhftc.m b/slhftc.m new file mode 100644 index 0000000..bb28e6c --- /dev/null +++ b/slhftc.m @@ -0,0 +1,55 @@ +function [qsen,qlat,tau,theta,A]=slhftc(ua,va,uo,vo,zr,Ta,zt,rh,zq,Pa,Ts) +% SLHFTC: computes sensible and latent heat flux following TOGA/COARE. +% [qsen,qlat,tau,theta,A]=SLHFTC(ua,va,uo,vo,zr,Ta,zt,rh,zq,ap,Ts) computes +% the sensible and latent heat fluxes into the ocean and the surface wind +% stress based on the Fairall et al (1996) COARE code. (see HFBULKTC for +% description). Assumes input series are either column or row vectors; +% zr, zt, and zq are fixed scalars, and rh and/or Pa may be scalars. +% The output variables are column vectors and the column matrix A. NOTE: +% user must decide if cool-skin and Webb corrections are to be included. +% +% INPUT: ua,va - east, north wind components [m/s] +% uo,vo - east, north ocean surface currents [m/s] +% zr - wind measurement height [m] +% Ta - air temperature [C] +% zt - air temperature measurement height [m] +% rh - relative humidity [%] +% zq - rh measurement height [m] +% Ts - sea surface temperature [C] +% Pa - air pressure [mb] +% +% OUTPUT: qsen - sensible heat flux [W/m^2] +% qlat - latent heat flux [W/m^2] +% tau - wind stress magnitude [Pa] +% theta - direction of wind stress [deg CCW from east] +% A - 12 column matrix of auxilary diagnostic outputs +% (see HFBULKTC for details) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 (contributed by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% compute relative velocity magnitude and direction +du=ua-uo; +dv=va-vo; +spd=sqrt(du.^2+dv.^2); +theta=(180./pi).*atan2(dv,du); + +% change column vector to row vector if necessary +[m n] = size(spd); +if m > n + spd = spd'; +end +% keep theta a column vector +[m n] = size(theta); +if n > m + theta = theta'; +end + +% compute fluxes +A = hfbulktc(spd,zr,Ta,zt,rh,zq,Pa,Ts); qsen=A(:,1); % no cool_skin, Webb correction +qlat=A(:,2); +tau=A(:,4); + + diff --git a/soradna1.m b/soradna1.m new file mode 100644 index 0000000..4c4ff61 --- /dev/null +++ b/soradna1.m @@ -0,0 +1,148 @@ +function [z,sorad]=soradna1(yd,yr,long,lat) +% SORADNA1: computes no-sky solar radiation and solar altitude. +% [z,sorad]=SORADNA1(yd,yr,long,lat) computes instantaneous values of +% solar radiation and solar altitude from yearday, year, and position +% data. It is put together from expressions taken from Appendix E in the +% 1978 edition of Almanac for Computers, Nautical Almanac Office, U.S. +% Naval Observatory. They are reduced accuracy expressions valid for the +% years 1800-2100. Solar declination computed from these expressions is +% accurate to at least 1'. The solar constant (1368.0 W/m^2) represents a +% mean of satellite measurements made over the last sunspot cycle (1979-1995) +% taken from Coffey et al (1995), Earth System Monitor, 6, 6-10. Assumes +% yd is either a column or row vector, the other input variables are scalars, +% OR yd is a scalar, the other inputs matrices. +% +% INPUT: yd - decimal yearday (e.g., 0000Z Jan 1 is 0.0) +% yr - year (e.g., 1995) +% long - longitude (west is positive!) [deg] +% lat - latitude [deg] +% +% OUTPUT: z - solar altitude [deg] +% sorad- no atmosphere solar radiation [W/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/28/98: version 1.1 (vectorized by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +% convert yd to column vector if necessary +[n,m]=size(yd); +if m > n + yd=yd'; +end + +% convert yearday to calender time +gtime=greg2(yd,yr); + +SC=gtime(:,6); +MN=fix(gtime(:,5)); +HR=fix(gtime(:,4)); +D=fix(gtime(:,3)); +M=fix(gtime(:,2)); +Y=fix(gtime(:,1)); + +% convert to new variables +LONG=long; +LAT=lat; + +% two options - either long/lat are vectors, time is a scalar + +if length(LONG)==1 & length(LAT)>1, + LONG=LONG(ones(size(LAT))); +elseif length(LONG)>1 & length(LAT)==1, + LAT=LAT(ones(size(LAT))); +end; + +if length(SC)==1, + osiz=ones(size(LONG)); + SC=SC(osiz); + MN=MN(osiz); + HR=HR(osiz); + D=D(osiz); + M=M(osiz); + Y=Y(osiz); +elseif length(LONG)==1, + LONG=LONG(ones(size(SC))); + LAT=LAT(ones(size(SC))); +end; + +% constants + DTR=3.14159265/180; + RTD=1./DTR; + +% compute Universal Time in hours + UT = HR+(MN+SC./60.)./60; + +% compute Julian ephemeris date in days (Day 1 is 1 Jan 4713 B.C.=-4712 Jan 1) + JED=367.*Y-fix(7.*(Y+fix((M+9)./12))./4)+fix(275.*M./9)+D+1721013 + UT./24; + +% compute interval in Julian centuries since 1900 + T=(JED-2415020.0)./36525; + +% compute mean anomaly of the sun + G=358.475833+35999.049750.*T-.000150.*T.^2; + NG=fix(G./360); + G=(G-NG.*360).*DTR; + +% compute mean longitude of sun + L=279.696678+36000.768920.*T+.000303.*T.^2; + NL=fix(L./360); + L=(L-NL.*360).*DTR; + +% compute mean anomaly of Jupiter + JUP=225.444651+2880.0.*T+154.906654.*T; + NJUP=fix(JUP/360); + JUP=(JUP-NJUP.*360).*DTR; + +% compute longitude of the ascending node of the moon's orbit + NM=259.183275-1800.*T-134.142008.*T+.002078.*T.^2; + NNM=fix(NM./360); + NM=(NM-NNM.*360+360).*DTR; + +% compute mean anomaly of Venus + V=212.603219+58320.*T+197.803875.*T+.001286.*T.^2; + NV=fix(V/360); + V=(V-NV.*360.).*DTR; + +% compute sun theta + THETA=.397930.*sin(L)+.009999.*sin(G-L)+.003334.*sin(G+L)... + -.000208.*T.*sin(L)+.000042.*sin(2.*G+L)-.000040.*cos(L)... + -.000039.*sin(NM-L)-.000030.*T.*sin(G-L)-.000014.*sin(2.*G-L)... + -.000010.*cos(G-L-JUP)-.000010.*T.*sin(G+L); + +% compute sun rho + RHO=1.000421-.033503.*cos(G)-.000140.*cos(2*G)... + +.000084.*T.*cos(G)-.000033.*sin(G-JUP)+.000027.*sin(2.*G-2.*V); + +% compute declination + DECL=asin(THETA./sqrt(RHO)); + +% compute equation of time (in seconds of time) (L in degrees) + L = 276.697+0.98564734.*(JED-2415020.0); + L = (L - 360.*fix(L./360.)).*DTR; + EQT = -97.8.*sin(L)-431.3.*cos(L)+596.6.*sin(2.*L)-1.9.*cos(2.*L)... + +4.0.*sin(3.*L)+19.3.*cos(3.*L)-12.7.*sin(4.*L); + EQT = EQT./60; + L = L.*RTD; + +% compute local hour angle + GHA = 15.*(UT-12.) + 15.*EQT./60; + LHA = GHA - LONG; + +% compute radius vector + RV=sqrt(RHO); + +% compute solar altitude + SZ=sin(DTR.*LAT).*sin(DECL)+cos(DTR.*LAT).*cos(DECL).*cos(DTR.*LHA); + z=RTD.*asin(SZ); + +% compute solar radiation outside atmosphere +[n,m]=size(z); +sorad=zeros(n,m); +ii=z>0; +sorad(ii)=(Solar_const./RV(ii).^2).*sin(DTR.*z(ii)); + diff --git a/spshftlp.m b/spshftlp.m new file mode 100644 index 0000000..f5e4467 --- /dev/null +++ b/spshftlp.m @@ -0,0 +1,30 @@ +function [sp2,ustar]=spshfttc(sp1,z1,z2) +% SPSHFTTC: adjusts wind speed from z1 to z2 following Large&Pond (1981). +% sp2 = SPSHFTLP(sp1,z1,z2) shifts the wind speed sp1 measured at z1 to +% z2 using the neutral drag coefficient given the wind speed and air +% temperature at height z following Large and Pond (1981), J. Phys. Oceanog., +% 11, 324-336. +% +% INPUT: sp1 - measured wind speed [m/s] +% z1 - measurement height [m] +% z2 - desired height [m] +% +% OUTPUT: sp2 - predicted wind speed [m/s] +% ustar - friction velocity [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/27/98: version 1.1 (revised to use CDNLP efficiently by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +% find cd and ustar +[cd10,sp10]=cdnlp(sp1,z1); + +ustar=sqrt(cd10).*sp10; + +sp2=sp10+ustar.*log(z2./10)/kappa; + diff --git a/spshfttc.m b/spshfttc.m new file mode 100644 index 0000000..7da205b --- /dev/null +++ b/spshfttc.m @@ -0,0 +1,34 @@ +function [sp2,ustar]=spshfttc(sp1,z1,z2,Ta) +% SPSHFTTC: adjusts wind speed from z1 to z2 following Smith (1988). +% sp2 = SPSHFTTC(sp1,z1,z2,Ta) shifts the wind speed sp1 measured at z1 to +% z2 using the neutral drag coefficient given the wind speed and air +% temperature at height z following Smith (1988), J. Geophys. Res., 93, +% 311-326. Assumes z1 and z2 scalars. Ta may be a constant. +% +% INPUT: sp1 - measured wind speed [m/s] +% z1 - measurement height [m] +% z2 - desired height [m] +% Ta - air temperature ([C] (optional) +% +% OUTPUT: sp2 - predicted wind speed [m/s] +% ustar - fiction velocity [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/27/98: version 1.1 (revised to use CDNTC efficiently by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% set constants +as_consts; + +if nargin==3, + Ta=Ta_default; +end; + +% find cd and ustar +[cd,sp10]=cdntc(sp1,z1,Ta); + +ustar=sqrt(cd).*sp10; + +sp2=sp10+ustar.*log(z2./10)/kappa; diff --git a/spshftve.m b/spshftve.m new file mode 100644 index 0000000..a8f7dd3 --- /dev/null +++ b/spshftve.m @@ -0,0 +1,30 @@ +function [sp2,ustar]=spshftve(sp1,z1,z2) +% SPSHFTVE: adjusts wind speed from z1 to z2 following Vera (1983). +% sp2 = SPSHFTVE(sp1,z1,z2) shifts the wind speed sp1 measured at z1 to +% z2 using the neutral drag law of Vera (1983) [see Large, Morzel, +% and Crawford (1995), J. Phys. Oceanog., 25, 2959-2971 (eqn. 8)]. +% Assumes z1 and z2 scalars. +% +% INPUT: sp1 - measured wind speed [m/s] +% z1 - measurement height [m] +% z2 - desired height of sp2 [m] +% +% OUTPUT: sp2 - predicted wind speed [m/s] +% ustar - friction velocity [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/27/98: version 1.1 (revised to use CDNVE efficiently by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% set constants +as_consts; + +% find cd and ustar +[cd10,sp10]=cdnve(sp1,z1); + +ustar=sqrt(cd10).*sp10; + +sp2=sp10+ustar.*log(z2./10)/kappa; + diff --git a/stresslp.m b/stresslp.m new file mode 100644 index 0000000..e8c3901 --- /dev/null +++ b/stresslp.m @@ -0,0 +1,30 @@ +function tau=stresslp(sp,z,rhoa) +% STRESSLP: computes neutral wind stress following Large and Pond (1981). +% tau = STRESSLP(sp,z,rhoa) computes the neutral wind stress given the wind +% speed at height z following Large and Pond (1981), J. Phys. Oceanog., +% 11, 324-336. Air density is an optional input, otherwise assumed to +% be constant (1.22 kg/m^3). +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% rhoa - air_density (optional) [kg/m^3] +% +% OUTPUT: tau - wind stress [N/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (revised by RP) +% 4/2/99: version 1.2 (optional air density added by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +if nargin == 2 + rhoa = rho_air; +end + +[cd,u10]=cdnlp(sp,z); +tau=rhoa.*(cd.*u10.^2); + diff --git a/stresstc.m b/stresstc.m new file mode 100644 index 0000000..13245a3 --- /dev/null +++ b/stresstc.m @@ -0,0 +1,34 @@ +function tau=stresstc(sp,z,Ta,rhoa) +% STRESSTC: computes the neutral wind stress following Smith (1988). +% tau = STRESSTC(sp,z,Ta,rhoa) computes the neutral wind stress given the +% wind speed and air temperature at height z following Smith (1988), +% J. Geophys. Res., 93, 311-326. Air temperature and density are optional +% inputs. +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% Ta - air temperature (optional) [C] +% rhoa - air density (optional) [kg/m^3] +% +% OUTPUT: tau - wind stress [N/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (revised by RP) +% 4/2/99: versin 1.2 (air density option added by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% load constants +as_consts; + +if nargin == 2, + Ta = Ta_default; + rhoa = rho_air; +elseif nargin == 3 + rhoa = rho_air; +end + +[cd,u10] = cdntc(sp,z,Ta); +tau = rhoa*(cd.*u10.^2); + diff --git a/stressve.m b/stressve.m new file mode 100644 index 0000000..a64ec5b --- /dev/null +++ b/stressve.m @@ -0,0 +1,32 @@ +function tau=stressve(sp,z,rhoa) +% STRESSVE: computes stress using Vera (1983) neutral drag law. +% tau = STRESSVE(sp,z,rhoa) computes the neutral wind stress given the wind +% speed at height z following Vera (1983) [see Large, Morzel, and Crawford +% (1995), J. Phys. Oceanog., 25, 2959-2971 (eqn. 8)]. Assumes z a fixed +% scalar. Air density is an optional input. +% +% INPUT: sp - wind speed [m/s] +% z - measurement height [m] +% rhoa - air density (optional) [kg/m^3] +% +% OUTPUT: tau - wind stress [N/m^2] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/26/98: version 1.1 (revised by RP) +% 4/2/99: version 1.2 (air density option added by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% load constants +as_consts; + +if nargin == 2 + rhoa = rho_air; +end + +[cd,u10]=cdnve(sp,z); +tau=rhoa*(cd.*u10.^2); + + + diff --git a/sunrise.m b/sunrise.m new file mode 100644 index 0000000..04d403c --- /dev/null +++ b/sunrise.m @@ -0,0 +1,61 @@ +function [rhr,rmin,shr,smin]=sunrise(mon,da,yr,lat,long) +% SUNRISE: computes sunrise and sunset times for specified day and location. +% [rhr,rmin,shr,smin] = SUNRISE(mon,da,yr,lat,long) computes the time +% of sunrise rhr:rmin and sunset shr:smin to the nearest minute in GMT +% for a calendar day(s) and a specified (scalar) position. +% +% INPUT: mon - month (e.g., Jan is 1) +% da - day (e.g., Jan 10th is 10) +% yr - year (e.g., 1995) +% lat - latitude [deg] +% long - longitude (west is positive) [deg] +% +% OUTPUT: rhr,rmin - sunrise in GMT hours and minutes +% shr,smin - sunset in GMT hours and minutes. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/28/98: version 1.1 (contributed by RP) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% convert calender time to julian yd +j=julianmd(yr,mon,da,0); +j0=julianmd(yr,1,1,0); +yd=j(:)-j0(:); + +% compute solar altitude for entire day +dt=1./2880; + +% we don't want abs(long)>180... +if long<-180, long=long+360; end; +if long>180, long=long-360; end; + +time=dt.*[0:2879]'+long/360; % have a whole day, beginning at midnight (near enough) +yday=yd(ones(1,2880),:)+time(:,ones(length(yd),1)); + +if length(yr)>1, + yr=yr(:,ones(1,2880))'; +end; + +[z,sorad]=soradna1(yday(:),yr(:),long,lat); + +z=reshape(z,2880,length(yd)); +sorad=reshape(sorad,2880,length(yd)); + +[ir,jr]=find(sorad(1:2879,:)==0 & sorad(2:2880,:)>0); +[is,js]=find(sorad(2:2880,:)==0 & sorad(1:2879,:)>0); + +srise=zeros(length(yd),1); +sset=any(sorad>0); + +srise(jr)=yday(ir+(jr-1)*2880); +sset(js) =yday(is+(js-1)*2880); + +rhr=fix(rem(srise,1)*24); +rmin=rem(rem(srise,1)*1440,60); +shr=fix(rem(sset,1)*24); +smin=rem(rem(sset,1)*1440,60); + + + + diff --git a/sw_tcond.m b/sw_tcond.m new file mode 100644 index 0000000..ae99496 --- /dev/null +++ b/sw_tcond.m @@ -0,0 +1,95 @@ +function tcond = sw_tcond(S,T,P) + +% SW_TCOND thermal conductivity +%=========================================================================== +% SW_TCOND $Revision: 0.0 $ $Date: 1998/01/19 $ +% Copyright (C) Ayal Anis 1998. +% +% USAGE: tcond = sw_tcond(S,T,P) +% +% DESCRIPTION: +% Calculates thermal conductivity of sea-water. +% Two options (Note: one option has to be remarked): +% 1) based on Caldwell's DSR 21:131-137 (1974) EQN. 9 +% 2) based on Catelli et al.'s DSR 21:311-3179(1974) EQN. 5 +% +% INPUT: (all must have same dimensions) +% S = salinity [psu (PSS-78) ] +% T = temperature [degree C (IPTS-68)] +% P = pressure [db] +% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) ) +% +% OUTPUT: +% tcond = thermal conductivity of sea-water [W/m/K] +% +% tcond(35,20,0)=0.5972 +% +% DISCLAIMER: +% This software is provided "as is" without warranty of any kind. +%========================================================================= + +% CALLER: general purpose +% CALLEE: none + +%---------------------- +% CHECK INPUT ARGUMENTS +%---------------------- +if nargin ~= 3 + error('sw_tcond.m: Must pass 3 parameters ') +end + +% CHECK S,T,P dimensions and verify consistent +[ms,ns] = size(S); +[mt,nt] = size(T); +[mp,np] = size(P); + + +% CHECK THAT S & T HAVE SAME SHAPE +if (ms~=mt) | (ns~=nt) + error('check_stp: S & T must have same dimensions') +end %if + +% CHECK OPTIONAL SHAPES FOR P +if mp==1 & np==1 % P is a scalar. Fill to size of S + P = P(1)*ones(ms,ns); +elseif np==ns & mp==1 % P is row vector with same cols as S + P = P( ones(1,ms), : ); % Copy down each column. +elseif mp==ms & np==1 % P is column vector + P = P( :, ones(1,ns) ); % Copy across each row +elseif mp==ms & np==ns % P is a matrix size(S) + % shape ok +else + error('check_stp: P has wrong dimensions') +end %if +[mp,np] = size(P); + + +% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN. +Transpose = 0; +if mp == 1 % row vector + P = P(:); + T = T(:); + S = S(:); + + Transpose = 1; +end + +%------ +% BEGIN +%------ + +% 1) Caldwell's option # 2 - simplified formula, accurate to 0.5% (eqn. 9) +tcond1 = 0.001365*(1+0.003*T-1.025e-5*T.^2+0.0653*(1e-4*P)-0.00029*S); % [cal/cm/C/sec] +tcond = tcond1*418.4; % [cal/cm/C/sec] ->[W/m/K] + +% 2) Castelli's option +%tcond2 = 100*(5.5286e-3+3.4025e-8*P+1.8364e-5*T-3.3058e-9*T.^3); % [W/m/K] +%tcond = tcond2 % [W/m/K] + +if Transpose + tcond = tcond'; +end + +return +%========================================================================= + diff --git a/sw_visc.m b/sw_visc.m new file mode 100644 index 0000000..af8eec2 --- /dev/null +++ b/sw_visc.m @@ -0,0 +1,71 @@ +function visc = sw_visc(S,T,P) + +% SW_VISC kinematic viscosity +%=========================================================================== +% SW_VISC $Revision: 0.0 $ $Date: 1998/01/19 $ +% Copyright (C) Ayal Anis 1998. +% +% USAGE: visc = sw_visc(S,T,P) +% +% DESCRIPTION: +% Calculates kinematic viscosity of sea-water. +% based on Dan Kelley's fit to Knauss's TABLE II-8 +% +% INPUT: (all must have same dimensions) +% S = salinity [psu (PSS-78) ] +% T = temperature [degree C (IPTS-68)] +% P = pressure [db] +% (P may have dims 1x1, mx1, 1xn or mxn for S(mxn) ) +% +% OUTPUT: +% visc = kinematic viscosity of sea-water [m^2/s] +% +% visc(40.,40.,1000.)=8.200167608E-7 +% +% DISCLAIMER: +% This software is provided "as is" without warranty of any kind. +%========================================================================= + +% CALLER: general purpose +% CALLEE: sw_dens.m + +%------------- +% CHECK INPUTS +%------------- +if nargin ~= 3 + error('sw_visc.m: Must pass 3 parameters ') +end + +% CHECK S,T dimensions and verify consistent +[ms,ns] = size(S); +[mt,nt] = size(T); + + % CHECK THAT S & T HAVE SAME SHAPE +if (ms~=mt) | (ns~=nt) + error('check_stp: S & T must have same dimensions') +end %if + +% LET sw_dens.m DO DIMENSION CHECKING FOR P + +% IF ALL ROW VECTORS ARE PASSED THEN LET US PRESERVE SHAPE ON RETURN. +Transpose = 0; +if ms == 1 % row vector + T = T(:); + S = S(:); + + Transpose = 1; +end %if + +%------ +% BEGIN +%------ + +visc = 1e-4*(17.91-0.5381*T+0.00694*T.^2+0.02305*S)./sw_dens(S,T,P); + +if Transpose + visc = visc'; +end %if + +return +%========================================================================= + diff --git a/swhf.m b/swhf.m new file mode 100644 index 0000000..c7eae0f --- /dev/null +++ b/swhf.m @@ -0,0 +1,42 @@ +function [qsw,alb]=swhf(yd,yr,long,lat,dsw) +% SWHF: computes net shortwave heat flux into the ocean and albedo. +% [qsw,alb]=SWHF(yd,yr,long,lat,dsw) computes the net shortwave heat +% flux into the ocean and albedo, using Payne (1972), J. Atm. Sci., 29, +% 959-970, to estimate the instantaneous albedo given the atmospheric +% transmittance (the ratio of measured insolation to the no-atmosphere +% insolation). +% +% INPUT: yd - decimal yearday (e.g., 0000Z Jan 1 is 0.0) +% yr - year (e.g., 1995) +% long - longitude [deg] +% lat - latitude [deg] +% dsw - (measured) insolation [W/m^2] +% +% OUTPUT: qsw - net shortwave heat flux [W/m^2] +% alb - albedo + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 7/29/99: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% compute sun altitude and no atm solar radiation +[sunalt, sorad]=soradna1(yd,yr,long,lat); + +% compute atm transmittance (note: trans=Inf when sorad=0) +trans=inf.*ones(size(sorad)); +j=find(sorad>0); +trans(j)=dsw(j)./sorad(j); + +% compute albedo (note: alb=NaN when trans>1 or sunalt<0) +alb=albedo(trans, sunalt); + +% compute net shortwave heat flux +% (note: qsw set equal to 0 when alb=NaN) +qsw=(1-alb).*dsw; +ind=find(isnan(qsw)); +qsw(ind)=zeros(size(ind)); + + + diff --git a/t_hfbulktc.m b/t_hfbulktc.m new file mode 100644 index 0000000..8cf9cfa --- /dev/null +++ b/t_hfbulktc.m @@ -0,0 +1,52 @@ +% T_HFBULKTC: a program to test hfbulktc.m using COARE test data +% T_HFBULKTC is a program to test hfbulktc.m using COARE test data +% (file test2_5b.dat). NOTE: Make sure that you use the COARE +% parameters located in the files as_consts.m and cool_skin.m when +% conducting the test. +% +% VARIBLES: +% +% 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)] +% dlw = downwelling (INTO water) longwave radiation [W/m^2] +% dsw = measured insolation [W/m^2] +% nsw = net shortwave radiation INTO the water [W/m^2] +% rain = rain rate [mm/hr] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 4/7/99: version 1.2 (contributed by AA) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% change directory to yours +load c:\flux_fsu\test\test2_5b.dat + +% parse data +ur = test2_5b(:,2); +zr = 15; +Ta = test2_5b(:,4); +zt = 15; +Pa = 1008*ones(size(ur)); +q = test2_5b(:,5); +rh = q./qsat(Ta,Pa)/10; +zq = 15; +Ts = test2_5b(:,14); +sal = 30*ones(size(ur)); +dsw = test2_5b(:,9); +nsw = 0.945*dsw; +dlw = test2_5b(:,10); +rain= test2_5b(:,11); + +i=[1:length(ur)]; +% i=[1:10]; +% NO cool-skin; compare to COARE output in file no_skin.out +A1 = hfbulktc(ur(i),zr,Ta(i),zt,rh(i),zq,Pa(i),Ts(i)); + +% YES cool-skin; compare to COARE output file yes_skin.out +A2 = hfbulktc(ur(i),zr,Ta(i),zt,rh(i),zq,Pa(i),Ts(i), ... + sal(i),dlw(i),dsw(i),nsw(i)); + diff --git a/test2_5b.dat b/test2_5b.dat new file mode 100644 index 0000000..b1692c5 --- /dev/null +++ b/test2_5b.dat @@ -0,0 +1,117 @@ +% datime GMT u tsea tair qair hsb hlb taub Rs Rl rain lat long MSP +921125132100 4.7 29 27.7 17.6 6 111 0.030 0 428 0 -1.727 156.065 29.15 +921125141200 4.1 29 27.7 17.7 5 96 0.023 0 429 0 -1.727 156.062 29.15 +921125150300 4.3 29 27.8 17.8 5 98 0.024 0 422 0 -1.729 156.057 29.151 +921125155500 4.7 29 27.8 17.6 5 110 0.029 0 412 0 -1.728 156.046 29.149 +921125164600 3.7 29 27.7 17.3 4 92 0.019 0 413 0 -1.728 156.035 29.161 +921125173800 5.1 29 27.3 17.6 9 123 0.036 0 419 0 -1.727 156.022 29.161 +921125182900 2.5 29 27.5 17.9 4 61 0.009 0 425 0 -1.725 156.014 29.161 +921125192000 3.8 29 27.9 17.7 4 89 0.019 50 427 0 -1.727 156.001 29.145 +921125201200 4.8 29 28 17.6 4 113 0.030 249 428 0 -1.728 155.995 29.136 +921125210500 3.9 29 27.8 18.3 5 84 0.021 386 422 0 -1.729 155.992 29.126 +921125221900 4.6 29.1 27.9 17.6 5 113 0.028 633 413 0 -1.728 155.987 29.143 +921125232700 5.2 29.2 28 18 6 124 0.036 881 409 0 -1.729 155.987 29.143 +921126001900 4.4 29.3 28.1 18.1 6 106 0.026 883 414 0 -1.731 155.983 29.199 +921126011000 4 29.3 28.1 18.6 5 90 0.022 853 420 0 -1.73 155.974 29.269 +921126022200 4.8 29.5 28.5 17.9 4 123 0.030 842 417 0 -1.728 155.969 29.258 +921126031300 4.8 29.5 28.4 18.1 5 119 0.030 779 411 0 -1.724 155.966 29.423 +921126041600 4.3 29.6 28.5 18.4 4 104 0.025 513 423 0 -1.721 155.96 29.586 +921126050700 4.4 29.5 28.6 18.9 3 94 0.025 362 422 0 -1.72 155.96 29.586 +921126055900 4.4 29.3 28.5 18.8 2 90 0.024 171 415 0 -1.721 155.959 29.5 +921126065500 4.1 29.3 28.5 17.9 2 96 0.022 26 410 0 -1.721 155.958 29.374 +921126074700 3.9 29.2 28.5 17.6 1 93 0.020 0 410 0 -1.721 155.954 29.374 +921126083800 3.2 29.3 28.5 17.8 1 78 0.014 0 410 0 -1.721 155.949 29.35 +921126093000 2.2 29.3 28.6 18 1 56 0.007 0 426 0 -1.721 155.945 29.387 +921126102100 1.9 29.3 28.5 18 1 50 0.006 0 414 0 -1.719 155.942 29.446 +921126111200 1.9 29.3 28.4 18 1 51 0.006 0 425 0 -1.717 155.939 29.418 +921126120800 2.5 29.2 28.4 18 1 60 0.009 0 417 0 -1.719 155.937 29.418 +921126130000 3 29.3 28.3 18.1 2 71 0.012 0 408 0 -1.719 155.938 29.352 +921126145300 2.7 29.2 28 18.7 3 58 0.011 0 412 0 -1.714 155.999 29.409 +921126154500 3 29.2 27.9 18.6 4 65 0.012 0 409 0 -1.715 156.024 29.309 +921126163600 3.1 29.2 27.9 18.5 4 68 0.013 0 413 0 -1.718 156.048 29.288 +921126172700 2.5 29.1 28 18.3 2 57 0.009 0 414 0 -1.721 156.068 29.274 +921126181900 1.9 29.1 28 18.2 2 47 0.006 0 409 0 -1.72 156.062 29.26 +921126191000 1.1 29.2 27.9 18.2 2 36 0.003 25 411 0 -1.719 156.055 29.274 +921126200100 0.8 29.2 28 18.4 1 30 0.002 110 414 0 -1.72 156.043 29.274 +921126213400 0.7 29.3 28.1 18.2 2 31 0.002 123 436 0 -1.722 156.017 29.293 +921126222600 2.1 29.2 27.3 18.4 5 55 0.007 116 445 0 -1.723 156.005 29.377 +921126231700 7.9 29.2 25.4 18.2 39 182 0.093 247 436 4.8 -1.724 156.001 29.31 +921127000900 6.6 29.3 25.9 17.7 30 175 0.064 635 418 0 -1.725 156.024 29.328 +921127010000 5.8 29.3 26.7 17.4 19 160 0.048 486 422 0 -1.725 156.056 29.33 +921127015100 4.8 29.3 26.9 17.5 14 129 0.032 356 429 0 -1.721 156.06 29.33 +921127024300 4.5 29.2 27 17.5 12 118 0.029 421 422 0 -1.721 156.048 29.361 +921127033400 5.2 29.2 26.8 17.9 15 129 0.037 191 434 0 -1.722 156.038 29.246 +921127042500 6.9 29.1 25.3 18.2 35 161 0.070 50 437 9.4 -1.723 156.036 29.225 +921127051700 5.7 29.1 25.1 17.5 32 154 0.048 72 434 1.6 -1.717 156.025 29.225 +921127060800 9.9 29.1 24.7 17.7 54 227 0.154 18 432 1.5 -1.718 156.012 29.237 +921127070000 5.2 29.1 25 17.6 29 137 0.040 6 429 0 -1.721 156.005 29.217 +921127075100 4.3 29.1 26.1 16.9 16 121 0.027 0 415 0 -1.716 155.992 29.217 +921127084200 3.7 29.1 26.7 16.3 10 112 0.020 0 414 0 -1.717 155.983 29.162 +921127093400 3.3 29.1 27.2 16.1 7 102 0.016 0 410 0 -1.716 155.977 29.154 +921127102500 3.9 29 27.5 16.2 6 113 0.021 0 410 0 -1.715 155.969 29.136 +921127113200 5 29 27.5 16.9 7 132 0.034 0 414 0 -1.716 155.96 29.136 +921127122600 5.9 29 27.6 17.3 8 147 0.048 0 415 0 -1.715 155.952 29.12 +921127134800 4.9 29 27.5 18.2 7 105 0.033 0 405 0 -1.717 155.946 29.108 +921127143900 5.2 29 27.6 18.3 7 110 0.036 0 406 0 -1.715 155.965 29.103 +921127160400 5.4 29 27.7 17.7 7 126 0.039 0 413 0 -1.717 156.008 29.125 +921127170000 4.4 29 27.7 17.7 5 102 0.026 0 416 0 -1.723 156.001 29.141 +921127175200 3.8 29 27.7 17.7 4 89 0.020 0 413 0 -1.73 155.995 29.124 +921127184300 3.1 29 27.7 17.6 4 75 0.013 2 406 0 -1.734 155.984 29.124 +921127193500 2.9 29 27.9 17.7 3 69 0.012 88 410 0 -1.72 155.998 29.123 +921127202600 2.1 29 27.9 17.8 2 54 0.007 276 410 0 -1.72 155.996 29.119 +921127211700 1.9 29.1 27.8 17.8 3 54 0.006 487 406 0 -1.724 155.988 29.119 +921127220900 1.6 29.2 27.7 17.7 3 52 0.005 670 407 0 -1.733 155.979 29.134 +921127230000 1.4 29.5 27.4 17.9 5 53 0.004 810 416 0 -1.738 155.98 29.138 +921127235200 1.1 30.1 27.4 17.8 7 56 0.003 917 410 0 -1.72 156.006 29.163 +921128004300 1.1 30.6 27.4 17.8 8 63 0.003 960 412 0 -1.718 155.999 29.136 +921128022100 1.5 30.8 27.7 17.9 9 74 0.005 939 413 0 -1.723 155.986 29.174 +921128042600 1.2 31 28.4 17.8 6 66 0.004 665 411 0 -1.716 156.001 29.271 +921128043600 1.3 31 28.7 17.5 5 69 0.004 493 410 0 -1.716 155.999 29.469 +921128052800 0.8 30.7 28.9 17.5 3 51 0.002 301 410 0 -1.72 155.993 29.234 +921128061900 0.6 30.6 28.9 17.6 2 44 0.002 125 411 0 -1.725 155.987 29.44 +921128075200 1.2 30.2 28.7 17.9 2 50 0.003 0 416 0 -1.717 155.999 29.273 +921128084300 1 30.1 28.7 17.9 2 45 0.003 0 417 0 -1.721 155.993 29.302 +921128093500 2.6 29.8 28.4 18.2 4 71 0.010 0 413 0 -1.721 155.991 29.302 +921128102600 1.8 29.5 28.5 17.5 2 56 0.005 0 420 0 -1.722 155.991 29.351 +921128111800 1 29.6 28.6 17.4 1 42 0.002 0 417 0 -1.721 155.991 29.242 +921128120900 2.8 29.3 28 18.8 3 61 0.011 0 411 0 -1.72 155.988 29.304 +921128130100 2.2 29.5 28 18.5 3 56 0.008 0 409 0 -1.719 155.984 29.559 +921128135200 1.6 29.4 28.1 18.1 2 47 0.004 0 407 0 -1.72 155.978 29.559 +921128144400 1.4 29.4 28.1 18 2 44 0.004 0 407 0 -1.721 155.969 29.288 +921128153500 1.9 29.2 28.1 17.8 2 52 0.006 0 407 0 -1.721 155.961 29.262 +921128162600 2.2 29.1 28 17.7 2 57 0.007 0 407 0 -1.721 155.952 29.262 +921128171800 2.3 29.1 28 18 2 56 0.008 0 406 0 -1.72 155.944 29.193 +921128180900 1.9 29 28 18 2 48 0.006 0 410 0 -1.72 155.936 29.134 +921128190100 2.3 29.1 28 17.8 2 58 0.008 14 407 0 -1.72 155.944 29.184 +921128195200 2.2 29.1 28.1 17.9 2 55 0.007 80 408 0 -1.717 155.968 29.169 +921128204400 1.6 29.2 28.1 17.5 2 49 0.005 109 414 0 -1.716 155.993 29.169 +921128213500 1.9 29.4 28.1 17.6 3 59 0.006 278 430 0 -1.713 156.017 29.399 +921128222600 1.2 29.6 27.2 18.5 6 46 0.003 614 425 0 -1.712 156.039 29.519 +921128231800 1.1 29.9 27.3 18.6 6 48 0.003 719 426 0 -1.713 156.049 29.519 +921129000900 0.5 30.4 27.1 18.4 7 47 0.002 930 415 0 -1.713 156.044 29.584 +921129010100 2.6 30.5 27.8 18.3 10 90 0.011 629 435 0 -1.708 156.038 29.492 +921129015200 2.5 30.1 27.7 17.9 9 84 0.011 838 418 0 -1.709 156.033 29.406 +921129024400 1.7 30.5 28.2 18 6 69 0.006 522 432 0 -1.711 156.025 29.406 +921129033500 2.1 30.3 28.5 17.7 5 77 0.008 745 413 0 -1.713 156.016 29.372 +921129042600 2.4 30.3 28.8 17.6 4 83 0.009 565 411 0 -1.714 156.005 29.23 +921129051800 2.1 30.1 28.8 17.8 3 70 0.007 209 430 0 -1.714 155.995 29.23 +921129060900 1.6 29.9 28.1 18.3 4 54 0.005 24 448 0 -1.717 155.982 29.315 +921129070100 3.7 30 26.6 17.9 17 112 0.021 1 445 6.5 -1.718 155.974 29.256 +921129075200 2.1 29.5 25.6 17.1 13 75 0.009 0 441 6.6 -1.717 155.966 29.459 +921129091700 2.2 29.8 27.3 17.1 7 76 0.008 0 410 0 -1.719 155.953 29.343 +921129100800 3.6 29.7 28.1 18.1 6 93 0.018 0 408 0 -1.719 155.947 29.331 +921129110000 3.2 29.4 28.2 18.4 3 74 0.014 0 409 0 -1.719 155.936 29.231 +921129115500 2.5 29.3 28.2 18.3 2 60 0.009 0 409 0 -1.718 155.945 29.231 +921129124700 2.6 29.2 28.1 18.3 2 60 0.010 0 407 0 -1.718 155.97 29.183 +921129133800 2.1 29.3 28 18.2 3 54 0.007 0 406 0 -1.717 155.995 29.402 +921129142900 2.1 29.5 28 18.2 3 57 0.007 0 405 0 -1.715 156.019 29.402 +921129153100 2.4 29.7 28 18.2 4 66 0.009 0 407 0 -1.715 156.046 29.627 +921129163800 3 29.6 28 18 5 79 0.013 0 411 0 -1.714 156.074 29.571 +921129173000 3.4 29.5 27.9 18.2 5 84 0.016 0 414 0 -1.716 156.064 29.494 +921129182100 3.1 29.4 28 18 4 77 0.014 0 417 0 -1.716 156.051 29.421 +921129191300 2.9 29.3 28 18 3 71 0.012 44 415 0 -1.715 156.038 29.421 +921129200400 3.2 29.3 27.8 18.1 5 78 0.014 213 418 0 -1.717 156.027 29.353 +921129205500 2.9 29.3 28.1 17.9 3 73 0.012 340 407 0 -1.717 156.018 29.352 +921129214700 2.4 29.4 28 17.9 4 67 0.009 597 405 0 -1.717 156.008 29.352 +921129223800 2.4 29.6 27.9 17.9 5 71 0.009 766 407 0 -1.717 155.999 29.369 +921129233000 2.4 29.8 27.8 17.8 7 77 0.009 900 411 0 -1.716 155.995 29.305 diff --git a/vapor.m b/vapor.m new file mode 100644 index 0000000..9c3919f --- /dev/null +++ b/vapor.m @@ -0,0 +1,27 @@ +function L=vapor(t) +% VAPOR calculates heat of evaporation for pure water +% L=VAPOR(t)computes the heat of evaporation for pure water. This can +% be used to compute the fresh water flux from latent heat flux. +% +% INPUT: t - water temperature [C] +% +% OUTPUT: L - heat of evaporation [J/kg] + +% Range of validity: 0 <= t <= 100 deg C. Check value: at t=100 deg C, +% L = 2.2566 x 10^6 J/kg. Reference: Landolt-Bornstein, Numerical Data +% and Functional Relationships in Science and Technology. New Series, +% Sundermann, J. (editor), vol. 3, subvol. a, Springer-Verlag, p. 256. +% No formulas are known to be available for the change of the heat of +% evaporation as function of salinity. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 1/22/99: version 1.0 (contributed by RO) +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +a0=2.50039e6; +a1=-2.3683e3; +a2=4.31e-1; +a3=-1.131e-2; +L=a0+a1*t+a2*t.*t+a3*t.*t.*t; + diff --git a/viscair.m b/viscair.m new file mode 100644 index 0000000..57a3ac4 --- /dev/null +++ b/viscair.m @@ -0,0 +1,17 @@ +function vis=viscair(Ta) +% VISCAIR: computes viscosity of air +% vis=VISCAIR(Ta) computes the kinematic viscosity of dry air as a +% function of air temperature following Andreas (1989), CRREL Report +% 89-11. +% +% INPUT: Ta - air temperature [C] +% +% OUTPUT: vis - air viscosity [m^2/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 3/8/97: version 1.0 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +vis = 1.326e-5*(1 + 6.542e-3*Ta + 8.301e-6*Ta.^2 - 4.84e-9*Ta.^3); + diff --git a/wavdist1.m b/wavdist1.m new file mode 100644 index 0000000..bebef0c --- /dev/null +++ b/wavdist1.m @@ -0,0 +1,49 @@ +function Ut=wavdist1(Ua,za,Hw) +% WAVDIST1: estimates wave effects on wind speed measured at za. +% Ut=WAVDIST1(Ua,za,Hw) computes the 'true' wind speed Ut at the +% measurement height za using the wind speed Ua measured at za and +% measured wave height Hw. +% +% INPUT: Ua - wind speed [m/s] +% za - wind measurement height [m] +% Hw - wave height [m] +% +% OUTPUT: Ut - 'true' wind speed [m/s] + +% WAVDIST1 computes Ut from Ua using the neutral log profile corrected +% for the effects of low-level distortion of the wind profile by surface +% waves following Large, Morzel, and Crawford (1995), J. Phys. Oceanog., +% 25, 2959-2971. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 5/5/97: version 1.0 +% 7/28/99: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +k=0.4; +z10=10; +A=log(z10./za)./k; + +% eliminate any Ua==0 +jj=find(Ua==0); +Ua(jj)=0.01.*ones(size(Ua(jj))); + +% compute uncorrected 10m wind speed +u10=Ua; % initial guess +for n=1:10; + ustar=sqrt(cdnve(u10).*u10.^2); + u10=Ua+ustar.*A; +end + +% compute corrected 10m wind speed +Ustar=ustar;U10=u10; % initial guesses +Za=za./Hw;Z10=z10./Hw; +for n=1:10; + Ustar=sqrt(cdnve(U10).*U10.^2); + U10=Ua+Ustar.*(log(z10./za)-omegalmc(Z10)+omegalmc(Za))./k; +end + +% compute 'true' wind speed at za using U10, Ustar +Ut=U10-Ustar.*A; + diff --git a/wavdist2.m b/wavdist2.m new file mode 100644 index 0000000..a6f7af1 --- /dev/null +++ b/wavdist2.m @@ -0,0 +1,35 @@ +function y=wavdist2(za) +% WAVDIST2: plots wave distortion effects on wind at za. +% WAVDIST2(za) plots the effects of wave distortion on the +% wind Ua measured at height za for the following significant +% wave heights Hw=[0:2:8] in m. +% +% INPUT: za - wind measurement height [m] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 5/5/97: version 1.0 +% 4/10/98: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Hw=[0:2:8]; +Ua=[0.01:.01:20]'; + +N=length(Hw); +M=length(Ua); +Ut=zeros(M,N); + +clg +hold on +for n=1:N + Ut=wavdist1(Ua,za,Hw(n)); + plot(Ua,Ut) +end + +title(['Predicted effects of wave distortion on wind speed at height ',num2str(za),' m']) +xlabel('Measured wind speed Ua (m/s)') +ylabel('Predicted wind speed Ut (m/s)') +text(10,2,'Wave height increment = 2 m') +grid + + diff --git a/wavedist.m b/wavedist.m new file mode 100644 index 0000000..38bdb66 --- /dev/null +++ b/wavedist.m @@ -0,0 +1,52 @@ +function [U10,delU]=wavedist(Ua,za,Hw) +% WAVEDIST: estimates wind speed distortion due to surface waves. +% [U10,delU]=WAVEDIST(Ua,za,Hw) computes the true 10m wind speed U10 +% using the wind speed Ua measured at the height za and measured wave +% height Hw and the neutral log profile corrected for the effects of +% low-level distortion of the wind profile by surface waves following +% Large, Morzel, and Crawford (1995), J. Phys. Ocean., 25, 2959-2971. +% +% INPUT: Ua - wind speed [m/s] +% za - wind measurement height [m] +% Hw - wave height [m] +% +% OUTPUT: U10 - true 10m wind speed [m/s] +% delU - difference between true and uncorrected +% 10m wind speed [m/s] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/31/98: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% get constants +as_consts; + +tol=.001; % change in u10 to stop iteration + +zs=10; % reference height + +Xia=za./Hw;Xis=zs./Hw; + +% compute uncorrected 10m wind speed and ustar (as initial guess in iteration) +[cd10,u10]=cdnve(Ua,za); +Ustar=sqrt(cd10).*u10; + +% compute corrected 10m wind speed +U10=u10; +U10o=0; +k=0; +while max(abs(U10-U10o))>tol & k<15, + U10o=U10; + k=k+1; + Ustar=sqrt(cdnve(U10,10).*U10.^2); + U10=Ua+Ustar.*(log(zs./za)-omegalmc(Xis)+omegalmc(Xia))./kappa; +end + +if k==15, + warning('Iteration may not have converged'); +end; + +delU=U10-u10; + + diff --git a/wdnotes.m b/wdnotes.m new file mode 100644 index 0000000..80a960c --- /dev/null +++ b/wdnotes.m @@ -0,0 +1,37 @@ +% WDNOTES: notes on estimating wave effects on wind measurements. +% The following set of mfiles can be used to correct the wind speed +% Ua measured at height za for the effects of the wave boundary layer +% following the empirical model presented by Large, Morzel, and +% Crawford (1995), J. Phys. Oceanog., 25, 2959-2971. In particular, +% an analytic expression was found for the omega function (omegalmc.m) +% shown in their Fig. 9b, which allows the 'true' wind speed (Ut10) +% and stress at 10m (assumed above the wave boundary layer height) +% to be computed using wavedist.m and the true wind speed (Uta) at the +% measurement height za using wavedis1.m. The Large et al model assumes +% neutral stability (reasonable for high winds and wave conditions) +% and uses a 10-m neutral drag law (cdnve.m) based on Vera (1983; +% unpublished manuscript). This drag law follows Large and Pond (1982) +% for winds above 10 m/s but increases at lower wind speeds like +% Smith (1987). The wave field is specified by the significant wave +% height Hw. +% +% To compute 'true' wind speed Uta at za given Hw, use +% Uta=wavedis1(Ua,za,Hw). +% +% To compute 'true' wind speed Ut at 10m given Hw, use +% [Ut10,(Ut10-U10)]=wavedist(Ua,za,Hw). +% +% To plot the predicted effects of wave distortion on the wind Ua +% measured at the height za for a range of significant wave heights +% Hw=[0:2:8] in m, use +% y=wavedis2(za). +% +% Subroutines called: +% y=omegalmc(x) +% cd10=cdnve(u10) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 7/28/99: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + diff --git a/yearday.m b/yearday.m new file mode 100644 index 0000000..4967ad5 --- /dev/null +++ b/yearday.m @@ -0,0 +1,47 @@ +function yd=yearday(mon,day,leapyr) +% YEARDAY: converts calender month and day into yearday. +% yd = YEARDAY(mon,day,leapyr) converts calender month and day into yearday. +% If year is not a leap year, then omit the third input variable. If year +% is a leap year, then enter leapyr=1. +% +% INPUT: mon - month +% day - day +% leapyr - set to 1 if it is a leap year, otherwise omit +% +% OUTPUT: yd - year day +% +% Examples: 15 March 1995 = yearday(3,15) = 74 +% 15 March 1996 = yearday(3,15,1) = 75 + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% 8/19/98: version 1.1 +% 8/5/99: version 2.0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +m=mon; +d=day; + +% compute yd for non-leap year +if m==1,y=d;end; +if m==2,y=31+d;end; +if m==3,y=59+d;end; +if m==4,y=90+d;end; +if m==5,y=120+d;end; +if m==6,y=151+d;end; +if m==7,y=181+d;end; +if m==8,y=212+d;end; +if m==9,y=243+d;end; +if m==10,y=273+d;end; +if m==11,y=304+d;end; +if m==12,y=334+d;end; +end +yd=y; +% adjust for leap year +if nargin > 2, + if (leapyr==1 & m>=3) + yd=yd+1; + end +end + + + diff --git a/yes_skin_out.txt b/yes_skin_out.txt new file mode 100644 index 0000000..0aadd97 --- /dev/null +++ b/yes_skin_out.txt @@ -0,0 +1,119 @@ +% Output from Fortran program bulk_v25b.f, with cool-skin turned ON and warm-layer +% turned OFF, and using input file test2_5b.dat with COARE data +% index xtime hsb hlb tub ts HF EF TAU T0 Webb RainF rain dt_cool dt_warm tk_pwp + 1 921125132100. 6. 111. .030 29.15 7. 116. .029 28.86 4. 0. .0 .29 .00 19.00 + 2 921125141200. 5. 96. .023 29.15 6. 100. .022 28.86 4. 0. .0 .29 .00 19.00 + 3 921125150300. 5. 98. .024 29.15 6. 102. .024 28.85 4. 0. .0 .30 .00 19.00 + 4 921125155500. 5. 110. .029 29.15 6. 114. .029 28.83 4. 0. .0 .31 .00 19.00 + 5 921125164600. 4. 92. .019 29.16 5. 96. .018 28.83 3. 0. .0 .33 .00 19.00 + 6 921125173800. 9. 123. .036 29.16 10. 128. .035 28.85 5. 0. .0 .31 .00 19.00 + 7 921125182900. 4. 61. .009 29.16 5. 64. .009 28.85 3. 0. .0 .31 .00 19.00 + 8 921125192000. 4. 89. .019 29.15 4. 92. .019 28.86 3. 0. .0 .28 .00 19.00 + 9 921125201200. 4. 113. .030 29.14 5. 117. .030 28.89 4. 0. .0 .25 .00 19.00 + 10 921125210500. 5. 84. .021 29.13 5. 87. .020 28.91 3. 0. .0 .21 .00 19.00 + 11 921125221900. 5. 113. .028 29.14 6. 114. .028 28.93 4. 0. .0 .21 .00 19.00 + 12 921125232700. 6. 124. .036 29.14 6. 122. .035 28.97 4. 0. .0 .17 .00 19.00 + 13 921126001900. 6. 106. .026 29.20 5. 104. .025 29.05 4. 0. .0 .15 .00 19.00 + 14 921126011000. 5. 90. .022 29.27 5. 90. .021 29.15 3. 0. .0 .12 .00 19.00 + 15 921126022200. 4. 123. .030 29.26 3. 115. .030 29.10 4. 0. .0 .16 .00 19.00 + 16 921126031300. 5. 119. .030 29.42 5. 117. .030 29.23 4. 0. .0 .19 .00 19.00 + 17 921126041600. 4. 104. .025 29.59 4. 104. .024 29.38 4. 0. .0 .21 .00 19.00 + 18 921126050700. 3. 94. .025 29.59 4. 97. .025 29.36 4. 0. .0 .22 .00 19.00 + 19 921126055900. 2. 90. .024 29.50 4. 95. .025 29.24 3. 0. .0 .26 .00 19.00 + 20 921126065500. 2. 96. .022 29.37 2. 98. .022 29.06 3. 0. .0 .32 .00 19.00 + 21 921126074700. 1. 93. .020 29.37 2. 98. .020 29.04 3. 0. .0 .33 .00 19.00 + 22 921126083800. 1. 78. .014 29.35 2. 79. .013 29.01 3. 0. .0 .34 .00 19.00 + 23 921126093000. 1. 56. .007 29.39 1. 58. .007 29.08 2. 0. .0 .31 .00 19.00 + 24 921126102100. 1. 50. .006 29.45 1. 53. .005 29.11 2. 0. .0 .34 .00 19.00 + 25 921126111200. 1. 51. .006 29.42 2. 53. .005 29.10 2. 0. .0 .31 .00 19.00 + 26 921126120800. 1. 60. .009 29.42 2. 64. .009 29.09 2. 0. .0 .33 .00 19.00 + 27 921126130000. 2. 71. .012 29.35 2. 72. .012 29.01 2. 0. .0 .34 .00 19.00 + 28 921126145300. 3. 58. .011 29.41 4. 62. .010 29.08 2. 0. .0 .33 .00 19.00 + 29 921126154500. 4. 65. .012 29.31 4. 67. .012 28.98 3. 0. .0 .33 .00 19.00 + 30 921126163600. 4. 68. .013 29.29 4. 70. .013 28.97 3. 0. .0 .32 .00 19.00 + 31 921126172700. 2. 57. .009 29.27 3. 60. .009 28.95 2. 0. .0 .33 .00 19.00 + 32 921126181900. 2. 47. .006 29.26 2. 50. .006 28.92 2. 0. .0 .34 .00 19.00 + 33 921126191000. 2. 36. .003 29.27 2. 37. .002 28.94 1. 0. .0 .33 .00 19.00 + 34 921126200100. 1. 30. .002 29.27 2. 31. .001 28.99 1. 0. .0 .28 .00 19.00 + 35 921126213400. 2. 31. .002 29.29 2. 32. .001 29.08 1. 0. .0 .22 .00 19.00 + 36 921126222600. 5. 55. .007 29.38 6. 57. .007 29.14 3. 0. .0 .23 .00 19.00 + 37 921126231700. 39. 182. .093 29.31 40. 187. .093 29.06 12. 26. 4.8 .25 .00 19.00 + 38 921127000900. 30. 175. .064 29.33 30. 176. .063 29.07 9. 0. .0 .25 .00 19.00 + 39 921127010000. 19. 160. .048 29.33 19. 161. .047 29.06 7. 0. .0 .27 .00 19.00 + 40 921127015100. 14. 129. .032 29.33 14. 130. .032 29.07 6. 0. .0 .26 .00 19.00 + 41 921127024300. 12. 118. .029 29.36 13. 123. .028 29.10 5. 0. .0 .26 .00 19.00 + 42 921127033400. 15. 129. .037 29.25 16. 131. .037 28.98 6. 0. .0 .26 .00 19.00 + 43 921127042500. 35. 161. .070 29.23 36. 166. .070 28.95 10. 50. 9.4 .27 .00 19.00 + 44 921127051700. 32. 154. .048 29.23 33. 158. .047 28.91 9. 9. 1.6 .31 .00 19.00 + 45 921127060800. 54. 227. .154 29.24 56. 234. .154 28.97 15. 9. 1.5 .26 .00 19.00 + 46 921127070000. 29. 137. .040 29.22 31. 141. .039 28.89 8. 0. .0 .33 .00 19.00 + 47 921127075100. 16. 121. .027 29.22 17. 124. .026 28.86 6. 0. .0 .36 .00 19.00 + 48 921127084200. 10. 112. .020 29.16 11. 113. .019 28.79 4. 0. .0 .37 .00 19.00 + 49 921127093400. 7. 102. .016 29.15 7. 103. .015 28.78 4. 0. .0 .38 .00 19.00 + 50 921127102500. 6. 113. .021 29.14 6. 116. .021 28.77 4. 0. .0 .37 .00 19.00 + 51 921127113200. 7. 132. .034 29.14 8. 136. .033 28.81 5. 0. .0 .33 .00 19.00 + 52 921127122600. 8. 147. .048 29.12 9. 152. .047 28.81 6. 0. .0 .31 .00 19.00 + 53 921127134800. 7. 105. .033 29.11 8. 109. .032 28.80 4. 0. .0 .31 .00 19.00 + 54 921127143900. 7. 110. .036 29.10 8. 113. .036 28.81 5. 0. .0 .30 .00 19.00 + 55 921127160400. 7. 126. .039 29.13 8. 130. .039 28.82 5. 0. .0 .30 .00 19.00 + 56 921127170000. 5. 102. .026 29.14 6. 106. .025 28.83 4. 0. .0 .31 .00 19.00 + 57 921127175200. 4. 89. .020 29.12 5. 92. .019 28.80 3. 0. .0 .32 .00 19.00 + 58 921127184300. 4. 75. .013 29.12 4. 77. .013 28.78 3. 0. .0 .35 .00 19.00 + 59 921127193500. 3. 69. .012 29.12 3. 72. .012 28.81 3. 0. .0 .32 .00 19.00 + 60 921127202600. 2. 54. .007 29.12 3. 56. .007 28.86 2. 0. .0 .26 .00 19.00 + 61 921127211700. 3. 54. .006 29.12 3. 54. .006 28.91 2. 0. .0 .21 .00 19.00 + 62 921127220900. 3. 52. .005 29.13 3. 51. .004 29.01 2. 0. .0 .13 .00 19.00 + 63 921127230000. 5. 53. .004 29.14 4. 49. .004 29.14 2. 0. .0 .00 .00 19.00 + 64 921127235200. 7. 56. .003 29.16 4. 44. .003 29.16 2. 0. .0 .00 .00 19.00 + 65 921128004300. 8. 63. .003 29.14 4. 44. .003 29.14 2. 0. .0 .00 .00 19.00 + 66 921128022100. 9. 74. .005 29.17 4. 50. .004 29.17 2. 0. .0 .00 .00 19.00 + 67 921128042600. 6. 66. .004 29.27 2. 43. .003 29.20 1. 0. .0 .07 .00 19.00 + 68 921128043600. 5. 69. .004 29.47 1. 47. .003 29.28 1. 0. .0 .18 .00 19.00 + 69 921128052800. 3. 51. .002 29.23 0. 33. .001 29.01 1. 0. .0 .22 .00 19.00 + 70 921128061900. 2. 44. .002 29.44 0. 31. .001 29.15 1. 0. .0 .29 .00 19.00 + 71 921128075200. 2. 50. .003 29.27 0. 38. .003 28.95 1. 0. .0 .32 .00 19.00 + 72 921128084300. 2. 45. .003 29.30 0. 35. .002 28.98 1. 0. .0 .32 .00 19.00 + 73 921128093500. 4. 71. .010 29.30 2. 62. .009 28.98 2. 0. .0 .33 .00 19.00 + 74 921128102600. 2. 56. .005 29.35 1. 54. .005 29.02 2. 0. .0 .33 .00 19.00 + 75 921128111800. 1. 42. .002 29.24 0. 38. .002 28.92 1. 0. .0 .32 .00 19.00 + 76 921128120900. 3. 61. .011 29.30 3. 61. .011 28.99 2. 0. .0 .32 .00 19.00 + 77 921128130100. 3. 56. .008 29.56 4. 57. .007 29.21 2. 0. .0 .35 .00 19.00 + 78 921128135200. 2. 47. .004 29.56 3. 49. .004 29.19 2. 0. .0 .37 .00 19.00 + 79 921128144400. 2. 44. .004 29.29 2. 43. .003 28.93 2. 0. .0 .35 .00 19.00 + 80 921128153500. 2. 52. .006 29.26 2. 53. .006 28.91 2. 0. .0 .36 .00 19.00 + 81 921128162600. 2. 57. .007 29.26 3. 60. .007 28.91 2. 0. .0 .36 .00 19.00 + 82 921128171800. 2. 56. .008 29.19 2. 58. .008 28.85 2. 0. .0 .35 .00 19.00 + 83 921128180900. 2. 48. .006 29.13 2. 50. .006 28.79 2. 0. .0 .34 .00 19.00 + 84 921128190100. 2. 58. .008 29.18 2. 59. .008 28.84 2. 0. .0 .35 .00 19.00 + 85 921128195200. 2. 55. .007 29.17 2. 56. .007 28.85 2. 0. .0 .32 .00 19.00 + 86 921128204400. 2. 49. .005 29.17 2. 49. .004 28.86 2. 0. .0 .31 .00 19.00 + 87 921128213500. 3. 59. .006 29.40 3. 59. .006 29.17 2. 0. .0 .23 .00 19.00 + 88 921128222600. 6. 46. .003 29.52 5. 45. .003 29.44 2. 0. .0 .08 .00 19.00 + 89 921128231800. 6. 48. .003 29.52 5. 44. .003 29.52 2. 0. .0 .00 .00 19.00 + 90 921129000900. 7. 47. .002 29.58 5. 36. .001 29.58 2. 0. .0 .00 .00 19.00 + 91 921129010100. 10. 90. .011 29.49 6. 70. .010 29.37 3. 0. .0 .12 .00 19.00 + 92 921129015200. 9. 84. .011 29.41 6. 72. .009 29.30 3. 0. .0 .10 .00 19.00 + 93 921129024400. 6. 69. .006 29.41 3. 53. .005 29.29 2. 0. .0 .12 .00 19.00 + 94 921129033500. 5. 77. .008 29.37 2. 62. .007 29.26 2. 0. .0 .11 .00 19.00 + 95 921129042600. 4. 83. .009 29.23 0. 64. .008 29.05 2. 0. .0 .18 .00 19.00 + 96 921129051800. 3. 70. .007 29.23 0. 55. .006 29.00 2. 0. .0 .23 .00 19.00 + 97 921129060900. 4. 54. .005 29.32 2. 46. .004 29.08 2. 0. .0 .24 .00 19.00 + 98 921129070100. 17. 112. .021 29.26 12. 95. .019 28.98 5. 32. 6.5 .28 .00 19.00 + 99 921129075200. 13. 75. .009 29.46 13. 75. .008 29.13 4. 41. 6.6 .33 .00 19.00 + 100 921129091700. 7. 76. .008 29.34 5. 69. .008 28.97 3. 0. .0 .37 .00 19.00 + 101 921129100800. 6. 93. .018 29.33 4. 85. .017 29.00 3. 0. .0 .33 .00 19.00 + 102 921129110000. 3. 74. .014 29.23 3. 71. .013 28.91 2. 0. .0 .32 .00 19.00 + 103 921129115500. 2. 60. .009 29.23 2. 59. .009 28.90 2. 0. .0 .33 .00 19.00 + 104 921129124700. 2. 60. .010 29.18 2. 60. .009 28.85 2. 0. .0 .33 .00 19.00 + 105 921129133800. 3. 54. .007 29.40 3. 56. .007 29.05 2. 0. .0 .36 .00 19.00 + 106 921129142900. 3. 57. .007 29.40 3. 56. .007 29.04 2. 0. .0 .36 .00 19.00 + 107 921129153100. 4. 66. .009 29.63 4. 65. .008 29.26 3. 0. .0 .36 .00 19.00 + 108 921129163800. 5. 79. .013 29.57 5. 79. .012 29.22 3. 0. .0 .35 .00 19.00 + 109 921129173000. 5. 84. .016 29.49 6. 84. .016 29.16 3. 0. .0 .33 .00 19.00 + 110 921129182100. 4. 77. .014 29.42 4. 78. .013 29.09 3. 0. .0 .33 .00 19.00 + 111 921129191300. 3. 71. .012 29.42 4. 74. .012 29.10 3. 0. .0 .32 .00 19.00 + 112 921129200400. 5. 78. .014 29.35 5. 79. .014 29.08 3. 0. .0 .27 .00 19.00 + 113 921129205500. 3. 73. .012 29.35 4. 74. .012 29.08 3. 0. .0 .27 .00 19.00 + 114 921129214700. 4. 67. .009 29.35 4. 66. .008 29.15 2. 0. .0 .20 .00 19.00 + 115 921129223800. 5. 71. .009 29.37 4. 68. .009 29.23 3. 0. .0 .14 .00 19.00 + 116 921129233000. 7. 77. .009 29.31 5. 69. .009 29.22 3. 0. .0 .09 .00 19.00