1 |
gforget |
1.1 |
function [uvcurl]=calc_UV_curl(u,v,putCurlOnTpoints, doMask); |
2 |
|
|
%[uvcurl]=calc_UV_curl(u,v,putCurlOnTpoints, [doMask]); |
3 |
|
|
%object: compute curl of a vector field |
4 |
|
|
%inputs: u,v is the vector field of interest |
5 |
|
|
% putCurlOnTpoints states wherther to return the curl |
6 |
|
|
% at curl points (0) or on tracer points (1) |
7 |
|
|
% doMask: 0 to computer curl at all points (default) or 1 to |
8 |
|
|
% mask so that curl is only computed at points that do not |
9 |
|
|
% touch land (set to zero where touches land) |
10 |
|
|
%output: [uvcurl] is the curl |
11 |
|
|
|
12 |
|
|
gcmfaces_global; |
13 |
|
|
|
14 |
|
|
if ~isfield(mygrid,'RAZfull'); |
15 |
|
|
%get full RAZ (incl. 'extra line and column') |
16 |
|
|
fprintf('\n calc_UV_curl requires loading (once) RAZfull \n'); |
17 |
|
|
fprintf('into mygrid, so we now call grid_load_native_RAZ \n\n'); |
18 |
|
|
grid_load_native_RAZ; |
19 |
|
|
end; |
20 |
|
|
|
21 |
|
|
if nargin < 4 |
22 |
|
|
doMask=0; |
23 |
|
|
end |
24 |
|
|
|
25 |
|
|
if doMask |
26 |
|
|
u(find(u==0))=NaN; v(find(v==0))=NaN; |
27 |
|
|
end |
28 |
|
|
|
29 |
|
|
%do the exchange: with sign changes for u and v, ... |
30 |
|
|
nn=1; |
31 |
|
|
[U,V]=exch_UV_N(u,v,nn); |
32 |
|
|
%... without sign changes for DXC and DYC |
33 |
|
|
[DXC,DYC]=exch_UV_N(mygrid.DXC,mygrid.DYC,nn); |
34 |
|
|
U=U.*abs(DXC); V=V.*abs(DYC); |
35 |
|
|
|
36 |
|
|
%compute the curl field |
37 |
|
|
uvcurl=NaN*mygrid.RAZfull; |
38 |
|
|
|
39 |
|
|
for iF=1:uvcurl.nFaces; |
40 |
|
|
%define ucur, vcur as u,v on current face, with points that dont |
41 |
|
|
%contribute to curl on that face removed. |
42 |
|
|
ucur=U{iF}(2:end,:); |
43 |
|
|
vcur=V{iF}(:,2:end); |
44 |
|
|
|
45 |
|
|
tmpcurl=ucur(:,1:end-1)-ucur(:,2:end); |
46 |
|
|
tmpcurl=tmpcurl-(vcur(1:end-1,:)-vcur(2:end,:)); |
47 |
|
|
|
48 |
|
|
%deal with corners where only 3 point contribute to curl |
49 |
|
|
if mod(iF,2)==1 %odd faces |
50 |
|
|
|
51 |
|
|
vc=-vcur(1,1); %- S-W corner |
52 |
|
|
tmpcurl(1,1) = (vcur(2,1)-ucur(1,2))+vc; |
53 |
|
|
|
54 |
|
|
vc=vcur(end,1); %- S-E corner |
55 |
|
|
tmpcurl(end,1) = (vc-ucur(end,2))-vcur(end-1,1); |
56 |
|
|
|
57 |
|
|
vc=vcur(end,end); %- N-E corner |
58 |
|
|
tmpcurl(end,end)=(vc-vcur(end-1,end))+ucur(end,end-1); |
59 |
|
|
|
60 |
|
|
vc=-vcur(1,end); %- N-W corner |
61 |
|
|
vc3=[vc vcur(2,end) ucur(1,end-1) vc vcur(2,end)]; |
62 |
|
|
n=(iF+1)/2; |
63 |
|
|
tmpcurl(1,end) = (vc3(n+2)+vc3(n+1))+vc3(n); |
64 |
|
|
|
65 |
|
|
else %even faces |
66 |
|
|
vc=-vcur(1,1); %- S-W corner |
67 |
|
|
tmpcurl(1,1) = (vcur(2,1)-ucur(1,2))+vc; |
68 |
|
|
|
69 |
|
|
vc=vcur(end,1); %- S-E corner |
70 |
|
|
n=iF/2; |
71 |
|
|
vc3=[-ucur(end,2) -vcur(end-1,1) vc -ucur(end,2) -vcur(end-1,1)]; |
72 |
|
|
tmpcurl(end,1) = (vc3(n)+vc3(n+1))+vc3(n+2); |
73 |
|
|
|
74 |
|
|
vc=vcur(end,end); %- N-E corner |
75 |
|
|
tmpcurl(end,end)=(ucur(end,end-1)+vc)-vcur(end-1,end); |
76 |
|
|
|
77 |
|
|
vc=-vcur(1,end); %- N-W corner |
78 |
|
|
tmpcurl(1,end) = (vcur(2,end)+vc)+ucur(1,end-1); |
79 |
|
|
end |
80 |
|
|
|
81 |
|
|
tmpcurl=tmpcurl./mygrid.RAZfull{iF}; |
82 |
|
|
|
83 |
|
|
%set points where NaN back to zero |
84 |
|
|
tmpcurl(isnan(tmpcurl))=0; |
85 |
|
|
|
86 |
|
|
%put to tracer points |
87 |
|
|
if putCurlOnTpoints; |
88 |
|
|
tmpcurl=1/4*(tmpcurl(1:end-1,2:end)+tmpcurl(1:end-1,1:end-1)+... |
89 |
|
|
tmpcurl(2:end,2:end)+tmpcurl(2:end,1:end-1)); |
90 |
|
|
end; |
91 |
|
|
uvcurl{iF}=tmpcurl; |
92 |
|
|
|
93 |
|
|
end; |
94 |
|
|
|