/[MITgcm]/MITgcm_contrib/gael/matlab_class/gcmfaces_smooth/diffsmooth2Drotated.m
ViewVC logotype

Diff of /MITgcm_contrib/gael/matlab_class/gcmfaces_smooth/diffsmooth2Drotated.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by gforget, Wed Feb 10 14:48:13 2010 UTC revision 1.2 by gforget, Wed Aug 1 00:27:13 2012 UTC
# Line 1  Line 1 
1  function [FLD]=diffsmooth2D(fld,dxLarge,dxSmall,fldRef);  function [FLD]=diffsmooth2Drotated(fld,dxLarge,dxSmall,fldRef);
2    %
3  %object: implementation (gforget@mit.edu) of a diffusive smoother (Weaver and Courtier 2001)  %object: slanted diffusive smoother (after Weaver and Courtier 2001)
4    %
5  %input: fld     field to be smoothed (masked with NaN)  %input: fld                  field to be smoothed (masked with NaN)
6  %       dxCsm,dyCsm     scale in first/second direction  %       dxLarge,dySmall  smoothing scale in direction of
7  %output:FLD     smoothed field  %                        weak,strong fldRef gradient
8    %       fldRef           tracer field which gradient defines
9  %asumption: dxCsm/dyCsm are assumed to be given at the positions of U/V points  %                        directions of strong,weak smoothing
10    %output:FLD                  smoothed field
11    %
12    %asumption: dxLarge/dxSmall are given at tracer points (not U/V points)
13    
14  global mygrid;  global mygrid;
15    
# Line 14  dxC=mygrid.DXC; dyC=mygrid.DYC; Line 17  dxC=mygrid.DXC; dyC=mygrid.DYC;
17  dxG=mygrid.DXG; dyG=mygrid.DYG;  dxG=mygrid.DXG; dyG=mygrid.DYG;
18  rA=mygrid.RAC;  rA=mygrid.RAC;
19    
 %compute the direction of main axis:  
 %%[dFLDdx,dFLDdy,ddFLDdxdx,ddFLDdydy]=calc_T_grad(fldRef,1);  
 %%dFLDn=sqrt(dFLDdx.^2+dFLDdy.^2);  
 %%ddFLDn=sqrt(ddFLDdxdx.^2+ddFLDdydy.^2);  
 [dFLDdx,dFLDdy]=calc_T_grad(fldRef,1);  
 dFLDn=sqrt(dFLDdx.^2+dFLDdy.^2);  
 if 0;  
    %FLDratio=dFLDn.*sqrt(mygrid.RAC)./max(fldRef,100);  
    FLDratio=dFLDn.*dxLarge./max(fldRef,100);  
 end;  
 cs=dFLDdy; sn=-dFLDdx;    
 %%cs=dFLDdx; sn=dFLDdy;  
 cs(dFLDn>0)=cs(dFLDn>0)./dFLDn(dFLDn>0);  
 sn(dFLDn>0)=sn(dFLDn>0)./dFLDn(dFLDn>0);  
 if 0;  
    FLDangle=fldRef;  
    for iF=1:fldRef.nFaces; FLDangle{iF}=atan2(sn{iF},cs{iF}); end;  
 end;  
   
20  %scale the diffusive operator:  %scale the diffusive operator:
21  tmp0=dxLarge./dxC; tmp0(isnan(fld))=NaN; tmp00=nanmax(tmp0);  tmp0=dxLarge./dxC; tmp0(isnan(fld))=NaN; tmp00=nanmax(tmp0);
22  tmp0=dxLarge./dyC; tmp0(isnan(fld))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);  tmp0=dxLarge./dyC; tmp0(isnan(fld))=NaN; tmp00=max([tmp00 nanmax(tmp0)]);
23  smooth2D_nbt=tmp00;  nbt=tmp00;
24  smooth2D_nbt=ceil(1.1*2*smooth2D_nbt^2);  nbt=ceil(1.1*2*nbt^2);
   
 smooth2D_dt=1;  
 smooth2D_T=smooth2D_nbt*smooth2D_dt;  
   
 smooth2D_kLarge=dxLarge.*dxLarge/smooth2D_T/2;  
 smooth2D_kSmall=dxSmall.*dxSmall/smooth2D_T/2;  
25    
26  smooth2D_Kux=cs.*cs.*smooth2D_kLarge+sn.*sn.*smooth2D_kSmall;  dt=1;
27  smooth2D_Kuy=cs.*sn.*(-smooth2D_kLarge+smooth2D_kSmall);  T=nbt*dt;
 smooth2D_Kvy=sn.*sn.*smooth2D_kLarge+cs.*cs.*smooth2D_kSmall;  
 smooth2D_Kvx=cs.*sn.*(-smooth2D_kLarge+smooth2D_kSmall);  
   
 %time-stepping loop:  
 FLD=fld;  
   
 for it=1:smooth2D_nbt;  
     if mod(it,ceil(smooth2D_nbt/50))==0; fprintf([num2str(it) '/' num2str(smooth2D_nbt) ' done\n']); end;  
   
     [dTdxAtU,dTdyAtV]=calc_T_grad(FLD,0);  
   
     dTdyAtU=dTdxAtU; dTdxAtV=dTdyAtV;  
     [dTdxAtU,dTdyAtV]=exch_UV_N(dTdxAtU,dTdyAtV);  
     for iF=1:FLD.nFaces;  
         msk=dTdxAtU{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;  
         tmp1=dTdyAtV{iF}; tmp1(isnan(tmp1))=0;  
         dTdyAtU{iF}=0.25*msk.*(tmp1(1:end-2,2:end-1)+tmp1(1:end-2,3:end)+...  
                 tmp1(2:end-1,2:end-1)+tmp1(2:end-1,3:end));  
   
         msk=dTdyAtV{iF}(2:end-1,2:end-1); msk(~isnan(msk))=1;  
         tmp1=dTdxAtU{iF}; tmp1(isnan(tmp1))=0;  
         dTdxAtV{iF}=0.25*msk.*(tmp1(2:end-1,1:end-2)+tmp1(3:end,1:end-2)+...  
                 tmp1(2:end-1,2:end-1)+tmp1(3:end,2:end-1));  
     end;  
     dTdxAtU=cut_T_N(dTdxAtU);  
     dTdyAtV=cut_T_N(dTdyAtV);  
   
     tmpU=dTdxAtU.*smooth2D_Kux+dTdyAtU.*smooth2D_Kuy;  
     tmpV=dTdyAtV.*smooth2D_Kvy+dTdxAtV.*smooth2D_Kvx;  
     [fldDIV]=calc_UV_div(tmpU,tmpV);  
     dFLDdt=smooth2D_dt*fldDIV./rA;  
     FLD=FLD-dFLDdt;  
28    
29  end;  %diffusion operator:
30    kLarge=dxLarge.*dxLarge/T/2;
31    kSmall=dxSmall.*dxSmall/T/2;
32    [Kux,Kuy,Kvx,Kvy]=diffrotated(kLarge,kSmall,fldRef);
33    
34    %setup problem:
35    myOp.dt=1;
36    myOp.nbt=nbt;
37    myOp.Kux=Kux;
38    myOp.Kuy=Kuy;
39    myOp.Kvx=Kvx;
40    myOp.Kvy=Kvy;
41    
42    %time step problem:
43    FLD=gcmfaces_timestep(myOp,fld);
44    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22