/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_UV_curl.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/gcmfaces_calc/calc_UV_curl.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Fri Sep 14 03:39:32 2012 UTC (12 years, 10 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65x, checkpoint65r, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, HEAD
- added gcmfaces_calc/calc_UV_curl.m : computes the curl of a vector field (from M.Buckley).

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    

  ViewVC Help
Powered by ViewVC 1.1.22