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	gebbie	1.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