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

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

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

Revision 1.2 - (hide annotations) (download)
Mon Aug 29 16:49:04 2011 UTC (13 years, 11 months ago) by gforget
Branch: MAIN
Changes since 1.1: +22 -65 lines 
>> clean up and speed up transport computations.

- added routines
gcmfaces_edge_mask		infer edge masks associated with an interior subdomain mask
gcmfaces_subset			extract subset of unmasked points from field (e.g. for transports)

- modified routines
gcmfaces_lines_zonal	use gcmfaces_edge_mask and make field names in
						LATS_MASKS more explicit (mskCint, mskCedge, etc.)
gcmfaces_lines_transp	same thing for section masks; also remove the
						longer arc masks (we only use the shorter arcs)
calc_zonmean_T			use new names (see gcmfaces_lines_zonal)
calc_zonmedian_T		use new names (see gcmfaces_lines_zonal)
calc_overturn			use new names (see gcmfaces_lines_zonal) and gcmfaces_subset
calc_MeridionalTransport use new names (see gcmfaces_lines_zonal) and gcmfaces_subset
calc_transports			use new names (see gcmfaces_lines_transp) and gcmfaces_subset
disp_transport			fix case of a single record
1 gforget 1.1 function [line_out]=gcmfaces_lines_transp(varargin);
2    
3     global mygrid;
4    
5     if nargin>0;
6     lonPairs=varargin{1};
7     latPairs=varargin{2};
8     names=varargin{3};
9     else;
10     lonPairs=[-68 -63];
11     latPairs=[-54 -66];
12     names={'Drake Passage'};
13     end;
14    
15     for iy=1:length(names);
16    
17     lonPair=lonPairs(iy,:);
18     latPair=latPairs(iy,:);
19     name=names{iy};
20    
21     %get carthesian coordinates:
22     %... of grid
23     lon=mygrid.XC; lat=mygrid.YC;
24     x=cos(lat*pi/180).*cos(lon*pi/180);
25     y=cos(lat*pi/180).*sin(lon*pi/180);
26     z=sin(lat*pi/180);
27     %... and of end points
28     x0=cos(latPair*pi/180).*cos(lonPair*pi/180);
29     y0=cos(latPair*pi/180).*sin(lonPair*pi/180);
30     z0=sin(latPair*pi/180);
31    
32     %get the rotation matrix:
33     %1) rotate around x axis to put first point at z=0
34     theta=atan2(-z0(1),y0(1));
35     R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]];
36     tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
37     x0=x1; y0=y1; z0=z1;
38     %2) rotate around z axis to put first point at y=0
39     theta=atan2(x0(1),y0(1));
40     R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]];
41     tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
42     x0=x1; y0=y1; z0=z1;
43     %3) rotate around y axis to put second point at z=0
44     theta=atan2(-z0(2),-x0(2));
45     R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]];
46     tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
47     x0=x1; y0=y1; z0=z1;
48    
49     %apply rotation to grid:
50     tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
51     tmp1=find(~isnan(tmpx));
52     tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1);
53     tmp2=[tmpx2';tmpy2';tmpz2'];
54     tmp3=R3*R2*R1*tmp2;
55     tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:);
56     tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2;
57     x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz);
58 gforget 1.2
59 gforget 1.1 %compute the great circle mask:
60 gforget 1.2 mskCint=1*(z>0);
61     [mskCedge,mskWedge,mskSedge]=gcmfaces_edge_mask(mskCint);
62    
63     %select the shorther segment:
64 gforget 1.1 for kk=1:3;
65     %select field to treat:
66     switch kk;
67 gforget 1.2 case 1; mm=mskCedge;
68     case 2; mm=mskWedge;
69     case 3; mm=mskSedge;
70 gforget 1.1 end;
71 gforget 1.2 %split in two segments:
72 gforget 1.1 theta=[];
73     theta(1)=atan2(y0(1),x0(1));
74     theta(2)=atan2(y0(2),x0(2));
75    
76     tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
77     tmptheta=atan2(tmpy,tmpx);
78 gforget 1.2 tmpm=convert2array(mm);
79 gforget 1.1 if theta(2)<0;
80     tmp00=find(tmptheta<=theta(2)); tmptheta(tmp00)=tmptheta(tmp00)+2*pi;
81     theta(2)=theta(2)+2*pi;
82     end;
83 gforget 1.2 %select the shorther segment:
84     if theta(2)-theta(1)<=pi;
85     tmpm(find(tmptheta>theta(2)|tmptheta<theta(1)))=NaN;
86     else;
87     tmpm(find(tmptheta<=theta(2)&tmptheta>=theta(1)))=NaN;
88     end;
89     mm=convert2array(tmpm);
90 gforget 1.1 %store result:
91     switch kk;
92 gforget 1.2 case 1; mskCedge=mm;
93     case 2; mskWedge=mm;
94     case 3; mskSedge=mm;
95 gforget 1.1 end;
96     end;
97    
98     %store so-defined line:
99 gforget 1.2 line_cur=struct('lonPair',lonPair,'latPair',latPair,'name',name,...
100     'mskCedge',mskCedge,'mskWedge',mskWedge,'mskSedge',mskSedge);
101 gforget 1.1
102     %add to lines:
103     if iy==1;
104     LINES_MASKS=line_cur;
105     else;
106     LINES_MASKS(iy)=line_cur;
107     end;
108    
109     end;
110    
111     mygrid.LINES_MASKS=LINES_MASKS;
112    
113    
114    
  ViewVC Help
Powered by ViewVC 1.1.22