| 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 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph