matlab_class/gcmfaces_calc/calc_T_grad.m

function [varargout]=calc_T_grad(fld,putGradOnTpoints);
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			global mygrid;
3
4			msk=fld; msk(~isnan(fld))=1;
5
6			FLD=exch_T_N(fld);
7			dFLDdx=fld; dFLDdy=fld;
8			for iF=1:FLD.nFaces;
9			tmpA=FLD{iF}(2:end-1,2:end-1);
10			tmpB=FLD{iF}(1:end-2,2:end-1);
11			dFLDdx{iF}=(tmpA-tmpB)./mygrid.DXC{iF};
12			tmpA=FLD{iF}(2:end-1,2:end-1);
13			tmpB=FLD{iF}(2:end-1,1:end-2);
14			dFLDdy{iF}=(tmpA-tmpB)./mygrid.DYC{iF};
15			end;
16
17			if nargout>2; %compute second order derivatives
18			[DX,DY]=exch_UV(dFLDdx,dFLDdy);
19			dFLDdxdx=fld; dFLDdydy=fld; dFLDdxdy=fld; dFLDdydx=fld;
20			for iF=1:FLD.nFaces;
21			tmpA=DX{iF}(2:end,:);
22			tmpB=DX{iF}(1:end-1,:);
23			dFLDdxdx{iF}=(tmpA-tmpB)./mygrid.DXF{iF};
24			tmpA=DY{iF}(:,2:end);
25			tmpB=DY{iF}(:,1:end-1);
26			dFLDdydy{iF}=(tmpA-tmpB)./mygrid.DYF{iF};
27			end;
28			end;
29
30			if putGradOnTpoints;
31			dFLDdx(isnan(dFLDdx))=0; dFLDdy(isnan(dFLDdy))=0;
32			dFLDdx0=dFLDdx; dFLDdy0=dFLDdy;
33			[dFLDdx,dFLDdy]=exch_UV(dFLDdx,dFLDdy);
34			for iF=1:FLD.nFaces;
35			dFLDdx{iF}=0.5*(dFLDdx{iF}(1:end-1,:)+dFLDdx{iF}(2:end,:));
36			dFLDdy{iF}=0.5*(dFLDdy{iF}(:,1:end-1)+dFLDdy{iF}(:,2:end));
37			end;
38			dFLDdx=dFLDdx.msk; dFLDdy=dFLDdy.msk;
39			%remakr: dFLDdxdx and dFLDdydy are on T points already
40			end;
41
42			varargout{1}=dFLDdx;
43			varargout{2}=dFLDdy;
44			if nargout>2;
45			varargout{3}=dFLDdxdx;
46			varargout{4}=dFLDdydy;
47			end;
48