matlab_class/gcmfaces_calc/calc_T_grad.m

function [varargout]=calc_T_grad(fld,putGradOnTpoints);
% CALC_T_GRAD(fld) computes horizontal gradients of a 2D field (fld)
%
% inputs:      fld is the two-dimensional field (at grid cell centers)
% [optional]   putGradOnTpoints states wherther to return the gradients
%                  at velocity points (0; default) or on tracer points (1)
% outputs:     [dFLDdx,dFLDdy] are the gradient fields
% [optional]   [...,dFLDdxdx,dFLDdydy] are the second order 
%                  derivatives at tracer points.
%

% development notes:
% - should be extended to 3D case

gcmfaces_global;

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.3	% CALC_T_GRAD(fld) computes horizontal gradients of a 2D field (fld)
3			%
4			% inputs: fld is the two-dimensional field (at grid cell centers)
5			% [optional] putGradOnTpoints states wherther to return the gradients
6			% at velocity points (0; default) or on tracer points (1)
7			% outputs: [dFLDdx,dFLDdy] are the gradient fields
8			% [optional] [...,dFLDdxdx,dFLDdydy] are the second order
9			% derivatives at tracer points.
10			%
11	gforget	1.2
12	gforget	1.3	% development notes:
13			% - should be extended to 3D case
14
15			gcmfaces_global;
16	gforget	1.1
17			msk=fld; msk(~isnan(fld))=1;
18
19			FLD=exch_T_N(fld);
20			dFLDdx=fld; dFLDdy=fld;
21			for iF=1:FLD.nFaces;
22			tmpA=FLD{iF}(2:end-1,2:end-1);
23			tmpB=FLD{iF}(1:end-2,2:end-1);
24			dFLDdx{iF}=(tmpA-tmpB)./mygrid.DXC{iF};
25			tmpA=FLD{iF}(2:end-1,2:end-1);
26			tmpB=FLD{iF}(2:end-1,1:end-2);
27			dFLDdy{iF}=(tmpA-tmpB)./mygrid.DYC{iF};
28			end;
29
30			if nargout>2; %compute second order derivatives
31			[DX,DY]=exch_UV(dFLDdx,dFLDdy);
32			dFLDdxdx=fld; dFLDdydy=fld; dFLDdxdy=fld; dFLDdydx=fld;
33			for iF=1:FLD.nFaces;
34			tmpA=DX{iF}(2:end,:);
35			tmpB=DX{iF}(1:end-1,:);
36			dFLDdxdx{iF}=(tmpA-tmpB)./mygrid.DXF{iF};
37			tmpA=DY{iF}(:,2:end);
38			tmpB=DY{iF}(:,1:end-1);
39			dFLDdydy{iF}=(tmpA-tmpB)./mygrid.DYF{iF};
40			end;
41			end;
42
43			if putGradOnTpoints;
44			dFLDdx(isnan(dFLDdx))=0; dFLDdy(isnan(dFLDdy))=0;
45			dFLDdx0=dFLDdx; dFLDdy0=dFLDdy;
46			[dFLDdx,dFLDdy]=exch_UV(dFLDdx,dFLDdy);
47			for iF=1:FLD.nFaces;
48			dFLDdx{iF}=0.5*(dFLDdx{iF}(1:end-1,:)+dFLDdx{iF}(2:end,:));
49			dFLDdy{iF}=0.5*(dFLDdy{iF}(:,1:end-1)+dFLDdy{iF}(:,2:end));
50			end;
51			dFLDdx=dFLDdx.msk; dFLDdy=dFLDdy.msk;
52			%remakr: dFLDdxdx and dFLDdydy are on T points already
53			end;
54
55			varargout{1}=dFLDdx;
56			varargout{2}=dFLDdy;
57			if nargout>2;
58			varargout{3}=dFLDdxdx;
59			varargout{4}=dFLDdydy;
60			end;
61