matlab_class/gcmfaces_calc/calc_T_grad.m

function [varargout]=calc_T_grad(fld,putGradOnTpoints);
%object:    compute horizontal gradients
%inputs:    fld is the 'tracer' field of interest
%           putGradOnTpoints states wherther to return the gradients
%               at velocity points (0) or on tracer points (1)
%output:    [dFLDdx,dFLDdy] are the gradient fields
%optional:  [...,dFLDdxdx,dFLDdydy] are the second order 
%               derivatives at tracer points.

global mygrid;

msk=fld; msk(~isnan(fld))=1;

FLD=exch_T_N(fld);
dFLDdx=fld; dFLDdy=fld;
for iF=1:FLD.nFaces;
   tmpA=FLD{iF}(2:end-1,2:end-1);
   tmpB=FLD{iF}(1:end-2,2:end-1);
   dFLDdx{iF}=(tmpA-tmpB)./mygrid.DXC{iF};
   tmpA=FLD{iF}(2:end-1,2:end-1);
   tmpB=FLD{iF}(2:end-1,1:end-2);
   dFLDdy{iF}=(tmpA-tmpB)./mygrid.DYC{iF};
end;

if nargout>2; %compute second order derivatives
   [DX,DY]=exch_UV(dFLDdx,dFLDdy);
   dFLDdxdx=fld; dFLDdydy=fld; dFLDdxdy=fld; dFLDdydx=fld;
   for iF=1:FLD.nFaces;
      tmpA=DX{iF}(2:end,:);
      tmpB=DX{iF}(1:end-1,:);
      dFLDdxdx{iF}=(tmpA-tmpB)./mygrid.DXF{iF};
      tmpA=DY{iF}(:,2:end);
      tmpB=DY{iF}(:,1:end-1);
      dFLDdydy{iF}=(tmpA-tmpB)./mygrid.DYF{iF};
   end;
end;

if putGradOnTpoints;
   dFLDdx(isnan(dFLDdx))=0; dFLDdy(isnan(dFLDdy))=0;
   dFLDdx0=dFLDdx; dFLDdy0=dFLDdy;
   [dFLDdx,dFLDdy]=exch_UV(dFLDdx,dFLDdy);
   for iF=1:FLD.nFaces;
      dFLDdx{iF}=0.5*(dFLDdx{iF}(1:end-1,:)+dFLDdx{iF}(2:end,:));
      dFLDdy{iF}=0.5*(dFLDdy{iF}(:,1:end-1)+dFLDdy{iF}(:,2:end));
   end;
   dFLDdx=dFLDdx.*msk; dFLDdy=dFLDdy.*msk;
%remakr: dFLDdxdx and dFLDdydy are on T points already
end;

varargout{1}=dFLDdx;
varargout{2}=dFLDdy;
if nargout>2;
   varargout{3}=dFLDdxdx;
   varargout{4}=dFLDdydy;
end;

1	gforget	1.1	function [varargout]=calc_T_grad(fld,putGradOnTpoints);
2	gforget	1.2	%object: compute horizontal gradients
3			%inputs: fld is the 'tracer' field of interest
4			% putGradOnTpoints states wherther to return the gradients
5			% at velocity points (0) or on tracer points (1)
6			%output: [dFLDdx,dFLDdy] are the gradient fields
7			%optional: [...,dFLDdxdx,dFLDdydy] are the second order
8			% derivatives at tracer points.
9
10	gforget	1.1	global mygrid;
11
12			msk=fld; msk(~isnan(fld))=1;
13
14			FLD=exch_T_N(fld);
15			dFLDdx=fld; dFLDdy=fld;
16			for iF=1:FLD.nFaces;
17			tmpA=FLD{iF}(2:end-1,2:end-1);
18			tmpB=FLD{iF}(1:end-2,2:end-1);
19			dFLDdx{iF}=(tmpA-tmpB)./mygrid.DXC{iF};
20			tmpA=FLD{iF}(2:end-1,2:end-1);
21			tmpB=FLD{iF}(2:end-1,1:end-2);
22			dFLDdy{iF}=(tmpA-tmpB)./mygrid.DYC{iF};
23			end;
24
25			if nargout>2; %compute second order derivatives
26			[DX,DY]=exch_UV(dFLDdx,dFLDdy);
27			dFLDdxdx=fld; dFLDdydy=fld; dFLDdxdy=fld; dFLDdydx=fld;
28			for iF=1:FLD.nFaces;
29			tmpA=DX{iF}(2:end,:);
30			tmpB=DX{iF}(1:end-1,:);
31			dFLDdxdx{iF}=(tmpA-tmpB)./mygrid.DXF{iF};
32			tmpA=DY{iF}(:,2:end);
33			tmpB=DY{iF}(:,1:end-1);
34			dFLDdydy{iF}=(tmpA-tmpB)./mygrid.DYF{iF};
35			end;
36			end;
37
38			if putGradOnTpoints;
39			dFLDdx(isnan(dFLDdx))=0; dFLDdy(isnan(dFLDdy))=0;
40			dFLDdx0=dFLDdx; dFLDdy0=dFLDdy;
41			[dFLDdx,dFLDdy]=exch_UV(dFLDdx,dFLDdy);
42			for iF=1:FLD.nFaces;
43			dFLDdx{iF}=0.5*(dFLDdx{iF}(1:end-1,:)+dFLDdx{iF}(2:end,:));
44			dFLDdy{iF}=0.5*(dFLDdy{iF}(:,1:end-1)+dFLDdy{iF}(:,2:end));
45			end;
46			dFLDdx=dFLDdx.msk; dFLDdy=dFLDdy.msk;
47			%remakr: dFLDdxdx and dFLDdydy are on T points already
48			end;
49
50			varargout{1}=dFLDdx;
51			varargout{2}=dFLDdy;
52			if nargout>2;
53			varargout{3}=dFLDdxdx;
54			varargout{4}=dFLDdydy;
55			end;
56