first commit

README.md

@@ -0,0 +1 @@
+# air-sea

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

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));
+
+

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                           

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
+
+
+

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

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;

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;