source: Roms_tools/Preprocessing_tools/smoothgrid.m @ 1

Last change on this file since 1 was 1, checked in by cholod, 11 years ago

import Roms_Agrif

File size: 5.1 KB
Line 
1function h = smoothgrid(h,maskr,hmin,hmax_coast,rmax,n_filter_deep_topo,n_filter_final)
2%
3%  Smooth the topography to get a maximum r factor = rmax
4%
5%  n_filter_deep_topo:
6%  Number of pass of a selective filter to reduce the isolated
7%  seamounts on the deep ocean.
8%
9%  n_filter_final:
10%  Number of pass of a single hanning filter at the end of the
11%  procedure to ensure that there is no 2DX noise in the
12%  topography.
13%
14%  Further Information: 
15%  http://www.brest.ird.fr/Roms_tools/
16
17%  This file is part of ROMSTOOLS
18%
19%  ROMSTOOLS is free software; you can redistribute it and/or modify
20%  it under the terms of the GNU General Public License as published
21%  by the Free Software Foundation; either version 2 of the License,
22%  or (at your option) any later version.
23%
24%  ROMSTOOLS is distributed in the hope that it will be useful, but
25%  WITHOUT ANY WARRANTY; without even the implied warranty of
26%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27%  GNU General Public License for more details.
28%
29%  You should have received a copy of the GNU General Public License
30%  along with this program; if not, write to the Free Software
31%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
32%  MA  02111-1307  USA
33%
34%  Copyright (c) 2005-2006 by Pierrick Penven
35%  e-mail:Pierrick.Penven@ird.fr 
36%
37%  Contributions of A. Shchepetkin (UCLA), P. Marchesiello (IRD)
38%                   and X. Capet (UCLA)
39%
40%  Updated    Aug-2006 by Pierrick Penven
41%
42%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43
44%
45% Cut the topography
46%
47h(h<hmin)=hmin;
48h(h>5000)=5000;
49%
50%
51% 1: Deep Ocean Filter
52%
53if n_filter_deep_topo>=1
54  disp(' Apply a filter on the Deep Ocean to remove the isolated seamounts :')
55  disp(['   ',num2str(n_filter_deep_topo),' pass of a selective filter.'])
56%
57%  Build a smoothing coefficient that is a linear function
58%  of a smooth topography.
59%
60  coef=h;
61  for i=1:8
62    coef=hanning_smoother(coef);             % coef is a smoothed bathy
63  end
64  coef=0.125*(coef./max(max(coef)));         % rescale the smoothed bathy
65
66  for i=1:n_filter_deep_topo;
67    h=hanning_smoother_coef2d(h,coef);       % smooth with avariable coef
68    h(maskr==0 & h>hmax_coast)=hmax_coast;
69  end
70end
71%
72%  Apply a selective filter on log(h) to reduce grad(h)/h.
73%
74disp(' Apply a selective filter on log(h) to reduce grad(h)/h :')
75h=rotfilter(h,maskr,hmax_coast,rmax);
76%
77%  Smooth the topography again to prevent 2D noise
78%
79if n_filter_final>1
80  disp(' Smooth the topography a last time to prevent 2DX noise:')
81  disp(['   ',num2str(n_filter_final),' pass of a hanning smoother.'])
82  for i=1:n_filter_final
83    h=hanning_smoother(h);
84    h(maskr==0 & h>hmax_coast)=hmax_coast;
85  end
86end
87%
88h(h<hmin)=hmin;
89%
90return
91%
92%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93function h=rotfilter(h,maskr,hmax_coast,rmax)
94%
95% Apply a selective filter on log(h) to reduce grad(h)/h.
96%
97[M,L]=size(h);
98Mm=M-1;
99Mmm=M-2;
100Lm=L-1;
101Lmm=L-2;
102cff=0.8;
103nu=3/16;
104[rx,ry]=rfact(h);
105r=max(max(max(rx)),max(max(ry)));
106h=log(h);
107hmax_coast=log(hmax_coast);
108i=0;
109while r>rmax
110  i=i+1;
111  cx=double(rx>cff*rmax);
112  cx=hanning_smoother(cx);
113  cy=double(ry>cff*rmax);
114  cy=hanning_smoother(cy);
115  fx=cx.*FX(h);
116  fy=cy.*FY(h);
117  h(2:Mm,2:Lm)=h(2:Mm,2:Lm)+nu*...
118             ((fx(2:Mm,2:Lm)-fx(2:Mm,1:Lmm))+...
119              (fy(2:Mm,2:Lm)-fy(1:Mmm,2:Lm)));
120  h(1,:)=h(2,:);
121  h(M,:)=h(Mm,:);
122  h(:,1)=h(:,2);
123  h(:,L)=h(:,Lm);
124  h(maskr==0 & h>hmax_coast)=hmax_coast;
125  [rx,ry]=rfact(exp(h));
126  r=max(max(max(rx)),max(max(ry)));
127end
128
129disp(['   ',num2str(i),' iterations - rmax = ',num2str(r)])
130h=exp(h);
131
132return
133%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134function [rx,ry]=rfact(h);
135[M,L]=size(h);
136Mm=M-1;
137Mmm=M-2;
138Lm=L-1;
139Lmm=L-2;
140rx=abs(h(1:M,2:L)-h(1:M,1:Lm))./(h(1:M,2:L)+h(1:M,1:Lm));
141ry=abs(h(2:M,1:L)-h(1:Mm,1:L))./(h(2:M,1:L)+h(1:Mm,1:L));
142
143return
144%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145function h=hanning_smoother(h);
146[M,L]=size(h);
147Mm=M-1;
148Mmm=M-2;
149Lm=L-1;
150Lmm=L-2;
151
152h(2:Mm,2:Lm)=0.125*(h(1:Mmm,2:Lm)+h(3:M,2:Lm)+...
153                       h(2:Mm,1:Lmm)+h(2:Mm,3:L)+...
154                       4*h(2:Mm,2:Lm));
155
156h(1,:)=h(2,:);
157h(M,:)=h(Mm,:);
158h(:,1)=h(:,2);
159h(:,L)=h(:,Lm);
160
161return
162%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163function h=hanning_smoother_coef2d(h,coef);
164[M,L]=size(h);
165Mm=M-1;
166Mmm=M-2;
167Lm=L-1;
168Lmm=L-2;
169h(2:Mm,2:Lm)=coef(2:Mm,2:Lm).*(h(1:Mmm,2:Lm)+h(3:M,2:Lm)+...
170                               h(2:Mm,1:Lmm)+h(2:Mm,3:L))...
171            +(1-4.*coef(2:Mm,2:Lm)).*h(2:Mm,2:Lm);
172h(1,:)=h(2,:);
173h(M,:)=h(Mm,:);
174h(:,1)=h(:,2);
175h(:,L)=h(:,Lm);
176
177return
178%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
179function fx=FX(h);
180[M,L]=size(h);
181Mm=M-1;
182Mmm=M-2;
183Lm=L-1;
184Lmm=L-2;
185
186fx(2:Mm,:)=(h(2:Mm,2:L)-h(2:Mm,1:Lm))*5/6 +...
187   (h(1:Mmm,2:L)-h(1:Mmm,1:Lm)+h(3:M,2:L)-h(3:M,1:Lm))/12;
188   
189fx(1,:)=fx(2,:);
190fx(M,:)=fx(Mm,:);
191
192return
193%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
194function fy=FY(h);
195[M,L]=size(h);
196Mm=M-1;
197Mmm=M-2;
198Lm=L-1;
199Lmm=L-2;
200
201fy(:,2:Lm)=(h(2:M,2:Lm)-h(1:Mm,2:Lm))*5/6 +...
202           (h(2:M,1:Lmm)-h(1:Mm,1:Lmm)+h(2:M,3:L)-h(1:Mm,3:L))/12;
203           
204fy(:,1)=fy(:,2);
205fy(:,L)=fy(:,Lm);
206
207return
208
Note: See TracBrowser for help on using the repository browser.