1 |
gforget |
1.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 |
|
|
|