matlab_class/gcmfaces_calc/gcmfaces_section.m

function [LO,LA,FLD,X,Y]=gcmfaces_section(lons,lats,fld,varargin);
%purpose: extract a great circle section (defined by two points) from a field
%         or a latitude circle (defined by one latitude)
%
%inputs:        lons/lats are the longitude/latitude vector
%               fld is the gcmfaces field (can incl. depth/time dimensions)
%optional:      sortByLon to sort point by longitude (default = 0 -> sort by latgitude) 
%outputs:       LO/LA is the vector of grid points longitude/latitude
%               FLD is the vector/matrix of grid point values (from fld)

gcmfaces_global;

if nargin>3; sortByLon=varargin{1}; else; sortByLon=0; end;

%check that LATS_MASKS has already been defined:
if ~isfield(mygrid,'LATS_MASKS');
    fprintf('one-time initialization of gcmfaces_lines_zonal: begin\n');
    gcmfaces_lines_zonal;
    fprintf('one-time initialization of gcmfaces_lines_zonal: end\n');
end;

if length(lats)==2;
    if lats(1)==-90&lats(2)==90; 
        error('need to specify latitude range under 180'); 
    end;
  line_cur=gcmfaces_lines_transp(lons,lats,{'tmp'});
elseif length(lats)==1;
  tmp1=abs(mygrid.LATS-lats);
  tmp2=find(tmp1==min(tmp1));
  tmp2=tmp2(1);
  line_cur=mygrid.LATS_MASKS(tmp2);
else;
  error('wrong specification of lons,lats');
end;
secP=find(line_cur.mskCedge==1);
secN=length(secP);

%lon/lat vectors:
LO=zeros(secN,1); LA=zeros(secN,1); 
%sections:
n3=max(size(fld{1},3),1); n4=max(size(fld{1},4),1); FLD=zeros(secN,n3,n4);
%counter:
ii0=0; 
for ff=1:secP.nFaces;
  tmp0=secP{ff}; [tmpI,tmpJ]=ind2sub(size(mygrid.XC{ff}),tmp0);
  tmp1=mygrid.XC{ff}; for ii=1:length(tmpI); LO(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
  tmp1=mygrid.YC{ff}; for ii=1:length(tmpI); LA(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
  tmp1=fld{ff}; for ii=1:length(tmpI); FLD(ii+ii0,:,:)=squeeze(tmp1(tmpI(ii),tmpJ(ii),:,:)); end;
  ii0=ii0+length(tmpI);
end;

%sort according to increasing latitude or longitude:
if sortByLon; 
  [tmp1,ii]=sort(LO); %sort according to increasing longitude
else;
  [tmp1,ii]=sort(LA); %sort according to increasing latitude
end;
LO=LO(ii); LA=LA(ii); FLD=FLD(ii,:,:);

%output axes for ploting with pcolor
nr=length(mygrid.RC);
nx=length(LO);
X=[]; Y=[];
if size(FLD,2)==nr&sortByLon; 
    X=LO'*ones(1,nr);
    Y=ones(nx,1)*mygrid.RC';
elseif size(FLD,2)==nr&~sortByLon;
    X=LA'*ones(1,nr);
    Y=ones(nx,1)*mygrid.RC';
elseif size(FLD,1)==nr&sortByLon; 
    X=LO';
    Y=ones(nx,1);
elseif size(FLD,2)==1&~sortByLon; 
    X=LA';
    Y=ones(nx,1);
end;
    
1	gforget	1.4	function [LO,LA,FLD,X,Y]=gcmfaces_section(lons,lats,fld,varargin);
2	gforget	1.1	%purpose: extract a great circle section (defined by two points) from a field
3	gforget	1.2	% or a latitude circle (defined by one latitude)
4	gforget	1.1	%
5			%inputs: lons/lats are the longitude/latitude vector
6			% fld is the gcmfaces field (can incl. depth/time dimensions)
7	gforget	1.2	%optional: sortByLon to sort point by longitude (default = 0 -> sort by latgitude)
8	gforget	1.4	%outputs: LO/LA is the vector of grid points longitude/latitude
9			% FLD is the vector/matrix of grid point values (from fld)
10	gforget	1.1
11	gforget	1.6	gcmfaces_global;
12
13	gforget	1.2	if nargin>3; sortByLon=varargin{1}; else; sortByLon=0; end;
14
15	gforget	1.6	%check that LATS_MASKS has already been defined:
16			if ~isfield(mygrid,'LATS_MASKS');
17			fprintf('one-time initialization of gcmfaces_lines_zonal: begin\n');
18			gcmfaces_lines_zonal;
19			fprintf('one-time initialization of gcmfaces_lines_zonal: end\n');
20			end;
21	gforget	1.1
22	gforget	1.2	if length(lats)==2;
23	gforget	1.5	if lats(1)==-90&lats(2)==90;
24			error('need to specify latitude range under 180');
25			end;
26	gforget	1.2	line_cur=gcmfaces_lines_transp(lons,lats,{'tmp'});
27			elseif length(lats)==1;
28			tmp1=abs(mygrid.LATS-lats);
29			tmp2=find(tmp1==min(tmp1));
30			tmp2=tmp2(1);
31			line_cur=mygrid.LATS_MASKS(tmp2);
32			else;
33			error('wrong specification of lons,lats');
34			end;
35			secP=find(line_cur.mskCedge==1);
36			secN=length(secP);
37	gforget	1.1
38			%lon/lat vectors:
39	gforget	1.4	LO=zeros(secN,1); LA=zeros(secN,1);
40	gforget	1.1	%sections:
41	gforget	1.4	n3=max(size(fld{1},3),1); n4=max(size(fld{1},4),1); FLD=zeros(secN,n3,n4);
42	gforget	1.1	%counter:
43			ii0=0;
44			for ff=1:secP.nFaces;
45			tmp0=secP{ff}; [tmpI,tmpJ]=ind2sub(size(mygrid.XC{ff}),tmp0);
46	gforget	1.4	tmp1=mygrid.XC{ff}; for ii=1:length(tmpI); LO(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
47			tmp1=mygrid.YC{ff}; for ii=1:length(tmpI); LA(ii+ii0)=tmp1(tmpI(ii),tmpJ(ii)); end;
48			tmp1=fld{ff}; for ii=1:length(tmpI); FLD(ii+ii0,:,:)=squeeze(tmp1(tmpI(ii),tmpJ(ii),:,:)); end;
49	gforget	1.1	ii0=ii0+length(tmpI);
50			end;
51
52	gforget	1.2	%sort according to increasing latitude or longitude:
53			if sortByLon;
54	gforget	1.4	[tmp1,ii]=sort(LO); %sort according to increasing longitude
55	gforget	1.2	else;
56	gforget	1.4	[tmp1,ii]=sort(LA); %sort according to increasing latitude
57	gforget	1.2	end;
58	gforget	1.4	LO=LO(ii); LA=LA(ii); FLD=FLD(ii,:,:);
59	gforget	1.1
60	gforget	1.4	%output axes for ploting with pcolor
61			nr=length(mygrid.RC);
62			nx=length(LO);
63			X=[]; Y=[];
64			if size(FLD,2)==nr&sortByLon;
65			X=LO'*ones(1,nr);
66			Y=ones(nx,1)*mygrid.RC';
67			elseif size(FLD,2)==nr&~sortByLon;
68			X=LA'*ones(1,nr);
69			Y=ones(nx,1)*mygrid.RC';
70			elseif size(FLD,1)==nr&sortByLon;
71			X=LO';
72			Y=ones(nx,1);
73			elseif size(FLD,2)==1&~sortByLon;
74			X=LA';
75			Y=ones(nx,1);
76			end;
77