MITgcm_contrib/gag_subduction_scripts/integrate_for_w.m

function [w] = integrate_for_w(u,v,dxg,dyg,raw, delZ )
%function [w] = integrate_for_w(u,v,dxg,dyg, raw, delZ )
%
% Get the vertical velocity from the horizontal velocity.
% Use the conservation of volume for this computation. 
% U and V are 3-dimensional.
% Following the MITgcm subroutine, integrate_for_w.F
%
% uncertain about the halo region.
%
% G. Gebbie, MIT-WHOI, 2003.
%

[nx ny nz] = size(u);

k=1:nz;
utrans(:,:,k) = u(:,:,k) .* dyg(:,:,ones(1,nz));
vtrans(:,:,k) = v(:,:,k) .* dxg(:,:,ones(1,nz));

%% is this the best way to overlap?
utrans(nx+1,:,k) = utrans(1,:,k);
vtrans(:,ny+1,k) = vtrans(:,1,k);

%w(:,:,23) = zeros(nx,ny,nz);

kbot = nz;
i=1:nx;
j=1:ny;
w(:,:,kbot) = - (utrans(i+1,j,kbot) - utrans(i,j,kbot) + ...
                 vtrans(i,j+1,kbot) - vtrans(i,j,kbot)) ...
                 .*(delZ(kbot)) ./raw(i,j);

for k=nz-1:-1:1
  w(:,:,k) = w(:,:,k+1)- ((utrans(i+1,j,k) - utrans(i,j,k) + ...
          vtrans(i,j+1,k) - vtrans(i,j,k)) .* (delZ(k) ./raw(i,j)));
end

return
1	function [w] = integrate_for_w(u,v,dxg,dyg,raw, delZ )
2	%function [w] = integrate_for_w(u,v,dxg,dyg, raw, delZ )
3	%
4	% Get the vertical velocity from the horizontal velocity.
5	% Use the conservation of volume for this computation.
6	% U and V are 3-dimensional.
7	% Following the MITgcm subroutine, integrate_for_w.F
8	%
9	% uncertain about the halo region.
10	%
11	% G. Gebbie, MIT-WHOI, 2003.
12	%
13
14	[nx ny nz] = size(u);
15
16	k=1:nz;
17	utrans(:,:,k) = u(:,:,k) .* dyg(:,:,ones(1,nz));
18	vtrans(:,:,k) = v(:,:,k) .* dxg(:,:,ones(1,nz));
19
20	%% is this the best way to overlap?
21	utrans(nx+1,:,k) = utrans(1,:,k);
22	vtrans(:,ny+1,k) = vtrans(:,1,k);
23
24	%w(:,:,23) = zeros(nx,ny,nz);
25
26	kbot = nz;
27	i=1:nx;
28	j=1:ny;
29	w(:,:,kbot) = - (utrans(i+1,j,kbot) - utrans(i,j,kbot) + ...
30	vtrans(i,j+1,kbot) - vtrans(i,j,kbot)) ...
31	.*(delZ(kbot)) ./raw(i,j);
32
33	for k=nz-1:-1:1
34	w(:,:,k) = w(:,:,k+1)- ((utrans(i+1,j,k) - utrans(i,j,k) + ...
35	vtrans(i,j+1,k) - vtrans(i,j,k)) .* (delZ(k) ./raw(i,j)));
36	end
37
38	return