matlab_class/gcmfaces_smooth/diffrotated.m

function [Kux,Kuy,Kvx,Kvy]=diffrotated(kLarge,kSmall,fldRef);
%
%object: compute slanted diffusive operator coefficients
%
%input: dxLarge,dySmall  smoothing scale in direction of 
%                        weak,strong fldRef gradient
%       fldRef           tracer field which gradient defines 
%                        directions of strong,weak smoothing
%output:Kux,Kuy,Kvx,Kvy  slanted diffusion operator coefficients
%
%asumption: dxLarge/dxSmall are given at tracer points (not U/V points)


gcmfaces_global;

%compute the direction of main axis:
%===================================

%1) gradient at cell center
[dTdxAtT,dTdyAtT]=calc_T_grad(fldRef,1);
[dTdxAtT,dTdyAtT]=exch_UV_N(dTdxAtT,dTdyAtT);
%2) diffusion coefficients at cell center
kLargeAtT=exch_T_N(kLarge); kSmallAtT=exch_T_N(kSmall);
%3) interpolate to U/V points
dTdxAtU=gcmfaces; dTdxAtV=gcmfaces;
dTdyAtU=gcmfaces; dTdyAtV=gcmfaces;
kLargeAtU=gcmfaces; kLargeAtV=gcmfaces;
kSmallAtU=gcmfaces; kSmallAtV=gcmfaces;
for iF=1:mygrid.nFaces;

    dTdxAtU{iF}=0.5*(dTdxAtT{iF}(1:end-2,2:end-1)+dTdxAtT{iF}(2:end-1,2:end-1));
    dTdyAtU{iF}=0.5*(dTdyAtT{iF}(1:end-2,2:end-1)+dTdyAtT{iF}(2:end-1,2:end-1));
    kLargeAtU{iF}=0.5*(kLargeAtT{iF}(1:end-2,2:end-1)+kLargeAtT{iF}(2:end-1,2:end-1));
    kSmallAtU{iF}=0.5*(kSmallAtT{iF}(1:end-2,2:end-1)+kSmallAtT{iF}(2:end-1,2:end-1));

    dTdxAtV{iF}=0.5*(dTdxAtT{iF}(2:end-1,1:end-2)+dTdxAtT{iF}(2:end-1,2:end-1));
    dTdyAtV{iF}=0.5*(dTdyAtT{iF}(2:end-1,1:end-2)+dTdyAtT{iF}(2:end-1,2:end-1));
    kLargeAtV{iF}=0.5*(kLargeAtT{iF}(2:end-1,1:end-2)+kLargeAtT{iF}(2:end-1,2:end-1));
    kSmallAtV{iF}=0.5*(kSmallAtT{iF}(2:end-1,1:end-2)+kSmallAtT{iF}(2:end-1,2:end-1));

end;


%compute diffusion operator : (rotated to the direction of main axis)
%===========================

%at U points
dFLDn=sqrt(dTdxAtU.^2+dTdyAtU.^2);
cs=dTdyAtU; sn=-dTdxAtU;
cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0); 
sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0);
if 1;
   FLDangleU=gcmfaces;
   for iF=1:mygrid.nFaces; FLDangleU{iF}=atan2(sn{iF},cs{iF}); end;
end;
Kux=cs.*cs.*kLargeAtU+sn.*sn.*kSmallAtU;
Kuy=cs.*sn.*(-kLargeAtU+kSmallAtU);
%at V points
dFLDn=sqrt(dTdxAtV.^2+dTdyAtV.^2);
cs=dTdyAtV; sn=-dTdxAtV;
cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0); 
sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0);
if 1;
   FLDangleV=fldRef;
   for iF=1:mygrid.nFaces; FLDangleV{iF}=atan2(sn{iF},cs{iF}); end;
end;
Kvy=sn.*sn.*kLargeAtV+cs.*cs.*kSmallAtV;
Kvx=cs.*sn.*(-kLargeAtV+kSmallAtV);

1	function [Kux,Kuy,Kvx,Kvy]=diffrotated(kLarge,kSmall,fldRef);
2	%
3	%object: compute slanted diffusive operator coefficients
4	%
5	%input: dxLarge,dySmall smoothing scale in direction of
6	% weak,strong fldRef gradient
7	% fldRef tracer field which gradient defines
8	% directions of strong,weak smoothing
9	%output:Kux,Kuy,Kvx,Kvy slanted diffusion operator coefficients
10	%
11	%asumption: dxLarge/dxSmall are given at tracer points (not U/V points)
12
13
14	gcmfaces_global;
15
16	%compute the direction of main axis:
17	%===================================
18
19	%1) gradient at cell center
20	[dTdxAtT,dTdyAtT]=calc_T_grad(fldRef,1);
21	[dTdxAtT,dTdyAtT]=exch_UV_N(dTdxAtT,dTdyAtT);
22	%2) diffusion coefficients at cell center
23	kLargeAtT=exch_T_N(kLarge); kSmallAtT=exch_T_N(kSmall);
24	%3) interpolate to U/V points
25	dTdxAtU=gcmfaces; dTdxAtV=gcmfaces;
26	dTdyAtU=gcmfaces; dTdyAtV=gcmfaces;
27	kLargeAtU=gcmfaces; kLargeAtV=gcmfaces;
28	kSmallAtU=gcmfaces; kSmallAtV=gcmfaces;
29	for iF=1:mygrid.nFaces;
30
31	dTdxAtU{iF}=0.5*(dTdxAtT{iF}(1:end-2,2:end-1)+dTdxAtT{iF}(2:end-1,2:end-1));
32	dTdyAtU{iF}=0.5*(dTdyAtT{iF}(1:end-2,2:end-1)+dTdyAtT{iF}(2:end-1,2:end-1));
33	kLargeAtU{iF}=0.5*(kLargeAtT{iF}(1:end-2,2:end-1)+kLargeAtT{iF}(2:end-1,2:end-1));
34	kSmallAtU{iF}=0.5*(kSmallAtT{iF}(1:end-2,2:end-1)+kSmallAtT{iF}(2:end-1,2:end-1));
35
36	dTdxAtV{iF}=0.5*(dTdxAtT{iF}(2:end-1,1:end-2)+dTdxAtT{iF}(2:end-1,2:end-1));
37	dTdyAtV{iF}=0.5*(dTdyAtT{iF}(2:end-1,1:end-2)+dTdyAtT{iF}(2:end-1,2:end-1));
38	kLargeAtV{iF}=0.5*(kLargeAtT{iF}(2:end-1,1:end-2)+kLargeAtT{iF}(2:end-1,2:end-1));
39	kSmallAtV{iF}=0.5*(kSmallAtT{iF}(2:end-1,1:end-2)+kSmallAtT{iF}(2:end-1,2:end-1));
40
41	end;
42
43
44	%compute diffusion operator : (rotated to the direction of main axis)
45	%===========================
46
47	%at U points
48	dFLDn=sqrt(dTdxAtU.^2+dTdyAtU.^2);
49	cs=dTdyAtU; sn=-dTdxAtU;
50	cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0);
51	sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0);
52	if 1;
53	FLDangleU=gcmfaces;
54	for iF=1:mygrid.nFaces; FLDangleU{iF}=atan2(sn{iF},cs{iF}); end;
55	end;
56	Kux=cs.cs.kLargeAtU+sn.sn.kSmallAtU;
57	Kuy=cs.sn.(-kLargeAtU+kSmallAtU);
58	%at V points
59	dFLDn=sqrt(dTdxAtV.^2+dTdyAtV.^2);
60	cs=dTdyAtV; sn=-dTdxAtV;
61	cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0);
62	sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0);
63	if 1;
64	FLDangleV=fldRef;
65	for iF=1:mygrid.nFaces; FLDangleV{iF}=atan2(sn{iF},cs{iF}); end;
66	end;
67	Kvy=sn.sn.kLargeAtV+cs.cs.kSmallAtV;
68	Kvx=cs.sn.(-kLargeAtV+kSmallAtV);
69