| 1 |
function [FLD]=diffsmooth2D(fld,dxLarge,dxSmall,fldRef); |
| 2 |
|
| 3 |
%object: implementation (gforget@mit.edu) of a diffusive smoother (Weaver and Courtier 2001) |
| 4 |
|
| 5 |
%input: fld field to be smoothed (masked with NaN) |
| 6 |
% dxCsm,dyCsm scale in first/second direction |
| 7 |
%output:FLD smoothed field |
| 8 |
|
| 9 |
%asumption: dxCsm/dyCsm are assumed to be given at the positions of U/V points |
| 10 |
|
| 11 |
global mygrid; |
| 12 |
|
| 13 |
dxC=mygrid.DXC; dyC=mygrid.DYC; |
| 14 |
dxG=mygrid.DXG; dyG=mygrid.DYG; |
| 15 |
rA=mygrid.RAC; |
| 16 |
|
| 17 |
%compute the direction of main axis: |
| 18 |
%%[dFLDdx,dFLDdy,ddFLDdxdx,ddFLDdydy]=calc_T_grad(fldRef,1); |
| 19 |
%%dFLDn=sqrt(dFLDdx.^2+dFLDdy.^2); |
| 20 |
%%ddFLDn=sqrt(ddFLDdxdx.^2+ddFLDdydy.^2); |
| 21 |
[dFLDdx,dFLDdy]=calc_T_grad(fldRef,1); |
| 22 |
dFLDn=sqrt(dFLDdx.^2+dFLDdy.^2); |
| 23 |
if 0; |
| 24 |
%FLDratio=dFLDn.*sqrt(mygrid.RAC)./max(fldRef,100); |
| 25 |
FLDratio=dFLDn.*dxLarge./max(fldRef,100); |
| 26 |
end; |
| 27 |
cs=dFLDdy; sn=-dFLDdx; |
| 28 |
%%cs=dFLDdx; sn=dFLDdy; |
| 29 |
cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0); |
| 30 |
sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0); |
| 31 |
if 0; |
| 32 |
FLDangle=fldRef; |
| 33 |
for iF=1:fldRef.nFaces; FLDangle{iF}=atan2(sn{iF},cs{iF}); end; |
| 34 |
end; |
| 35 |
|
| 36 |
%scale the diffusive operator: |
| 37 |
tmp0=dxLarge./dxC; tmp0(isnan(fld))=NaN; tmp00=nanmax(tmp0); |
| 38 |
tmp0=dxLarge./dyC; tmp0(isnan(fld))=NaN; tmp00=max([tmp00 nanmax(tmp0)]); |
| 39 |
smooth2D_nbt=tmp00; |
| 40 |
smooth2D_nbt=ceil(1.1*2*smooth2D_nbt^2); |
| 41 |
|
| 42 |
smooth2D_dt=1; |
| 43 |
smooth2D_T=smooth2D_nbt*smooth2D_dt; |
| 44 |
|
| 45 |
smooth2D_kLarge=dxLarge.*dxLarge/smooth2D_T/2; |
| 46 |
smooth2D_kSmall=dxSmall.*dxSmall/smooth2D_T/2; |
| 47 |
|
| 48 |
smooth2D_Kux=cs.*cs.*smooth2D_kLarge+sn.*sn.*smooth2D_kSmall; |
| 49 |
smooth2D_Kuy=cs.*sn.*(-smooth2D_kLarge+smooth2D_kSmall); |
| 50 |
smooth2D_Kvy=sn.*sn.*smooth2D_kLarge+cs.*cs.*smooth2D_kSmall; |
| 51 |
smooth2D_Kvx=cs.*sn.*(-smooth2D_kLarge+smooth2D_kSmall); |
| 52 |
|
| 53 |
%time-stepping loop: |
| 54 |
FLD=fld; |
| 55 |
|
| 56 |
for it=1:smooth2D_nbt; |
| 57 |
if mod(it,ceil(smooth2D_nbt/50))==0; fprintf([num2str(it) '/' num2str(smooth2D_nbt) ' done\n']); end; |
| 58 |
|
| 59 |
[dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0); |
| 60 |
|
| 61 |
dTdyAtU=dTdxAtU; dTdxAtV=dTdyAtV; |
| 62 |
[dTdxAtU,dTdyAtV]=exch_UV_N(dTdxAtU,dTdyAtV); |
| 63 |
for iF=1:FLD.nFaces; |
| 64 |
msk=dTdxAtU{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1; |
| 65 |
tmp1=dTdyAtV{iF}; tmp1(isnan(tmp1))=0; |
| 66 |
dTdyAtU{iF}=0.25*msk.*(tmp1(1:end-2,2:end-1)+tmp1(1:end-2,3:end)+... |
| 67 |
tmp1(2:end-1,2:end-1)+tmp1(2:end-1,3:end)); |
| 68 |
|
| 69 |
msk=dTdyAtV{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1; |
| 70 |
tmp1=dTdxAtU{iF}; tmp1(isnan(tmp1))=0; |
| 71 |
dTdxAtV{iF}=0.25*msk.*(tmp1(2:end-1,1:end-2)+tmp1(3:end,1:end-2)+... |
| 72 |
tmp1(2:end-1,2:end-1)+tmp1(3:end,2:end-1)); |
| 73 |
end; |
| 74 |
dTdxAtU=cut_T_N(dTdxAtU); |
| 75 |
dTdyAtV=cut_T_N(dTdyAtV); |
| 76 |
|
| 77 |
tmpU=dTdxAtU.*smooth2D_Kux+dTdyAtU.*smooth2D_Kuy; |
| 78 |
tmpV=dTdyAtV.*smooth2D_Kvy+dTdxAtV.*smooth2D_Kvx; |
| 79 |
[fldDIV]=calc_UV_div(tmpU,tmpV); |
| 80 |
dFLDdt=smooth2D_dt*fldDIV./rA; |
| 81 |
FLD=FLD-dFLDdt; |
| 82 |
|
| 83 |
end; |
| 84 |
|
| 85 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph