📄 cbepm.m
字号:
function [BEp, BEm]=cBEpm(g, f, sigma, nu, kappa, z, h)
%
%Evaluate the theoretical vertical profiles (or Bottom Ekman spiral profiles)
%of tidal currents in the two rotary directions driven by half-unit of sea
%surface gradients in the two directions respectively. Eddy viscosity is
%assumed as vertically invariant. See tidal_ellipse.ps for more details.
%
%
%inputs:
%
% g, the gravity acceleration,
% f, the Coriolis parameter,
% nu, the eddy viscosity
% kappa, the bottom frictional coefficient
% z, the vertical coordinates, can be a vector but must be
% within [0 -h];
% h, the water depth, must be positive.
%
% Note: except for z, all other inputs must be scalars.
%
%outputs:
%
% BEp and BEm, the same dimensions of z, the outputs for the vertical
% velocity profiles driven respectively by a unit of sea
% surface slope in the positive rotation direction and negative
% rotation direction for when the eddy viscosity is vertically
% invariant. See the associated document for more details.
if length(g)>1 | length(f)>1 | length(sigma)>1 | ...
length(nu)>1 | length(kappa)>1 | length(h)>1
error('inputs of g,f,sigma, nu, kappa, and h should be all scalars!')
end
if any(z/h>0) | any(z/h<-1)
disp('z must be negative and must be within [0 -h]')
end
delta_e = sqrt(2*nu/f); %Ekman depth
alpha = (1+i)/delta_e*sqrt(1+sigma/f);
beta = (1+i)/delta_e*sqrt(1-sigma/f);
BEp = get_BE(g, alpha, h, z, nu, kappa);
BEm = get_BE(g, beta, h, z, nu, kappa);
%subfunction
%______________________________________________
function BE=get_BE(g, alpha, h, z, nu, kappa)
z = z(:);
z_h = z/h;
ah = alpha*h;
az = alpha*z;
ah2 = ah*2;
anu_k = alpha*nu/kappa;
nu_kh = nu/(kappa*h);
if abs(ah) < 1 %series solution
T = 10;
C = -g*h*h/(nu*(1+anu_k*tanh_v5_2(ah)))*2;
A1 = (1-z_h.*z_h)/2+nu_kh;
B1 = exp(-ah)/(1+exp(-ah2));
B = B1;
series_sum=A1*B1;
for t = 2:T;
t2=2*t;
A = (1-z_h.^t2)./t2+nu_kh;
B = B*ah*ah/(t2-1)/(t2-2);
series_sum = series_sum+A*B;
end
BE = C*series_sum;
else %finite solution
c = -g*h*h/nu;
denom=(exp(az-ah)+exp(-(az+ah)))./(1+exp(-2*ah));
% =cosh(az)/cosh(ah);
%but this a better way to evaluate it.
numer=1+anu_k*tanh_v5_2(ah);
% BE=c*(1-denom/numer);
BE = c*((1-denom/numer)/(ah*ah));
end
%end of subfunction
%
%Note tanh_v5_2 is a copy of tanh from Matlab v5.2, which has worked well!
%It seems that Matlab v5.3 has some bug(s) in tanh function! It cannot deal
%with large argument. try z=7.7249e02*(1+i), tanh(z) and tanh_v5_2(z) to
%see the difference.
%Authorship Copyright:
%
% The author of this program retains the copyright of this program, while
% you are welcome to use and distribute this program as long as you credit
% the author properly and respect the program name itself. Particularly,
% you are expected to retain the original author's name in this original
% version of the program or any of its modified version that you might make.
% You are also expected not to essentially change the name of the programs
% except for adding possible extension for your own version you might create,
% e.g. ap2ep_xx is acceptable. Any suggestions are welcome and enjoying my
% program(s)!
%
%
%Author Info:
%_______________________________________________________________________
% Zhigang Xu, Ph.D.
% (pronounced as Tsi Gahng Hsu)
% Research Scientist
% Coastal Circulation
% Bedford Institute of Oceanography
% 1 Challenge Dr.
% P.O. Box 1006 Phone (902) 426-2307 (o)
% Dartmouth, Nova Scotia Fax (902) 426-7827
% CANADA B2Y 4A2 email zhigangx@emerald.bio.dfo.ca
% zhigang_xu_98@yahoo.com
%_______________________________________________________________________
%
%Release Date: Nov. 2000
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -