| 1 |
function [vort,z6t]=calc_vort_cs(u3d,v3d,dxC,dyC,rAz); |
| 2 |
% [vort,z6t]=calc_vort_cs(u3d,v3d,dxC,dyC,rAz); |
| 3 |
% compute vorticity (3rd component) on CS-grid. |
| 4 |
% assume u3d(nc*6,nc,*), v3d(nc*6,nc,*), dxC(nc*6,nc), dyC(nc*6,nc) |
| 5 |
% and deal with: rAz(nc*6,nc) or rAz(nc*6*nc+2) |
| 6 |
% output is provided with 2 shapes: |
| 7 |
% vort(nc*6*nc+2,*) = compressed form; |
| 8 |
% z6t(nc+1,nc+1,*,6) = face splitted. |
| 9 |
% $Header: /u/gcmpack/MITgcm/utils/cs_grid/calc_vort_cs.m,v 1.1 2005/09/06 17:38:06 jmc Exp $ |
| 10 |
% $Name: $ |
| 11 |
|
| 12 |
dims=size(u3d); nx=dims(1); nc=dims(2); ncp=nc+1; n2p=nc+2; nPg=nx*nc ; |
| 13 |
if nx == 6*nc, flag=1; else flag=0; end |
| 14 |
if size(v3d) ~= dims, flag=0; end |
| 15 |
if size(dxC) ~= dims(1:2), flag=0; end |
| 16 |
if size(dyC) ~= dims(1:2), flag=0; end |
| 17 |
if flag == 0, |
| 18 |
fprintf(' Error in size of input arrays\n'); |
| 19 |
vort=0; return |
| 20 |
end |
| 21 |
if length(dims) == 2, nr=1; else nr=prod(dims(3:end)); end |
| 22 |
|
| 23 |
siZ=prod(size(rAz)); |
| 24 |
if siZ == nPg+2, |
| 25 |
rAz=reshape(rAz,[nPg+2 1]); |
| 26 |
aZ=split_Z_cub(rAz); |
| 27 |
elseif siZ == nPg, |
| 28 |
rAz=reshape(rAz,[nPg 1]); |
| 29 |
aZc=zeros(nPg+2,1); aZc(1:nPg,:)=rAz; aZc(nPg+1)=rAz(1); aZc(nPg+2)=rAz(1); |
| 30 |
aZ=split_Z_cub(aZc); |
| 31 |
else |
| 32 |
fprintf(' Error in size of rAz input array\n'); |
| 33 |
vort=0; return |
| 34 |
end |
| 35 |
|
| 36 |
u3d=reshape(u3d,[nx nc nr]); |
| 37 |
v3d=reshape(v3d,[nx nc nr]); |
| 38 |
|
| 39 |
[u6t,v6t]=split_UV_cub(u3d,v3d,1,2); |
| 40 |
[d6x,d6y]=split_UV_cub(dxC,dyC,0,2); |
| 41 |
if nr > 1, |
| 42 |
u6t=permute(u6t,[1 2 4 3]); |
| 43 |
v6t=permute(v6t,[1 2 4 3]); |
| 44 |
end |
| 45 |
z6t=zeros(ncp,ncp,6,nr); |
| 46 |
|
| 47 |
for k=1:nr, |
| 48 |
vv1=d6x.*u6t(:,:,:,k); |
| 49 |
vv2=d6y.*v6t(:,:,:,k); |
| 50 |
z6t(:,:,:,k)=vv2([2:n2p],:,:)-vv2([1:ncp],:,:); |
| 51 |
z6t(:,:,:,k)=z6t(:,:,:,k)-(vv1(:,[2:n2p],:)-vv1(:,[1:ncp],:)); |
| 52 |
%- corner: the quick way: |
| 53 |
% z6t(1,1,:,k) = vv2(2,1,:) -vv2(1,1,:) -vv1(1,2,:); |
| 54 |
% z6t(1,ncp,:,k) = vv2(2,ncp,:)-vv2(1,ncp,:)+vv1(1,ncp,:); |
| 55 |
% z6t(ncp,1,:,k) = vv2(n2p,1,:)-vv2(ncp,1,:)-vv1(ncp,2,:); |
| 56 |
% z6t(ncp,ncp,:,k)=vv2(n2p,ncp,:)-vv2(ncp,ncp,:)+vv1(ncp,ncp,:); |
| 57 |
%- corner: add the 3 terms always in the same order |
| 58 |
% to get the same truncation on the 3 faces |
| 59 |
for n=1:3, |
| 60 |
f=2*n-1; %- odd face number |
| 61 |
vc=-vv2(1,1,f); %- S-W corner |
| 62 |
z6t(1,1,f,k) = (vv2(2,1,f)+vc) -vv1(1,2,f); |
| 63 |
vc=+vv2(n2p,1,f); %- S-E corner |
| 64 |
z6t(ncp,1,f,k) =-(vv2(ncp,1,f)+vv1(ncp,2,f))+vc; |
| 65 |
vc=+vv2(n2p,ncp,f); %- N-E corner |
| 66 |
z6t(ncp,ncp,f,k)=(vv1(ncp,ncp,f)+vc)-vv2(ncp,ncp,f); |
| 67 |
vc=-vv2(1,ncp,f); %- N-W corner |
| 68 |
vc3=[vv1(1,ncp,f) vc vv2(2,ncp,f) vv1(1,ncp,f) vc]; |
| 69 |
z6t(1,ncp,f,k) = (vc3(n+2)+vc3(n+1))+vc3(n); |
| 70 |
f=2*n; %- even face number |
| 71 |
vc=-vv2(1,1,f); %- S-W corner |
| 72 |
z6t(1,1,f,k) = (vc -vv1(1,2,f))+vv2(2,1,f); |
| 73 |
vc=+vv2(n2p,1,f); %- S-E corner |
| 74 |
vc3=[-vv1(ncp,2,f) -vv2(ncp,1,f) vc -vv1(ncp,2,f) -vv2(ncp,1,f)]; |
| 75 |
z6t(ncp,1,f,k) = (vc3(n)+vc3(n+1))+vc3(n+2); |
| 76 |
vc=+vv2(n2p,ncp,f); %- N-E corner |
| 77 |
z6t(ncp,ncp,f,k)=(vc-vv2(ncp,ncp,f))+vv1(ncp,ncp,f); |
| 78 |
vc=-vv2(1,ncp,f); %- N-W corner |
| 79 |
z6t(1,ncp,f,k) = (vv1(1,ncp,f)+vv2(2,ncp,f))+vc; |
| 80 |
end |
| 81 |
%- divide by rAz: |
| 82 |
z6t(:,:,:,k)=z6t(:,:,:,k)./aZ; |
| 83 |
end |
| 84 |
|
| 85 |
%- put in compressed form: |
| 86 |
vort=zeros(nPg+2,nr); |
| 87 |
% extract the interior |
| 88 |
%vv3=z6t(1:nc,1:nc,:,k); |
| 89 |
%vv3=permute(vv3,[1 3 2]); |
| 90 |
%vort([1:nPg],k)=reshape(vv3,[nPg 1]); |
| 91 |
vort([1:nPg],:)=reshape(permute(z6t(1:nc,1:nc,:,:),[1 3 2 4]),[nPg nr]); |
| 92 |
% put back the 2 missing corners (N.W corner of 1rst face & S.E corner of 2nd face) |
| 93 |
vort(nPg+1,:)=z6t(1,ncp,1,:); |
| 94 |
vort(nPg+2,:)=z6t(ncp,1,2,:); |
| 95 |
|
| 96 |
%- back into final shape: |
| 97 |
z6t=permute(z6t,[1 2 4 3]); |
| 98 |
if length(dims) == 2, |
| 99 |
vort=squeeze(vort); |
| 100 |
z6t=squeeze(z6t); |
| 101 |
else |
| 102 |
vort=reshape(vort,[nPg+2 dims(3:end)]); |
| 103 |
z6t=reshape(z6t,[ncp ncp dims(3:end) 6]); |
| 104 |
end |
| 105 |
|
| 106 |
%----------------- |
| 107 |
return |
Parent Directory
| 
Revision Log
| 
Revision Graph