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	function [varargout]=calc_T_grad(fld,putGradOnTpoints);
2	% 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
12	% development notes:
13	% - should be extended to 3D case
14
15	gcmfaces_global;
16
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