1 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
2 | % |
---|
3 | % Build a ROMS grid file |
---|
4 | % |
---|
5 | % |
---|
6 | % Pierrick Penven, IRD, 2002. |
---|
7 | % |
---|
8 | % Version of 10-Oct-2002 |
---|
9 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
10 | clear all |
---|
11 | close all |
---|
12 | %%%%%%%%%%%%%%%%%%%%% USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%%% |
---|
13 | % |
---|
14 | % Title |
---|
15 | % |
---|
16 | title='TEST'; |
---|
17 | config='test'; |
---|
18 | % |
---|
19 | % Grid file name |
---|
20 | % |
---|
21 | grdname='test_grd.nc'; |
---|
22 | % |
---|
23 | % Slope parameter (r=grad(h)/h) maximum value for topography smoothing |
---|
24 | % |
---|
25 | rtarget=0.15; |
---|
26 | % |
---|
27 | % Grid dimensions: |
---|
28 | % lonmin : Minimum longitude [degree east] |
---|
29 | % lonmax : Maximum longitude [degree east] |
---|
30 | % latmin : Minimum latitude [degree north] |
---|
31 | % latmax : Maximum latitude [degree north] |
---|
32 | if config=='test' |
---|
33 | % |
---|
34 | % Test |
---|
35 | lon1=43.9; |
---|
36 | lat1=8.7; |
---|
37 | lon2=34.1; |
---|
38 | lat2=31.2; |
---|
39 | lon3=46.9; |
---|
40 | lat3=10.7; |
---|
41 | |
---|
42 | end |
---|
43 | % |
---|
44 | % Grid resolution [degree] |
---|
45 | % |
---|
46 | dl=1/4; |
---|
47 | % |
---|
48 | % Minimum depth [m] |
---|
49 | % |
---|
50 | hmin=10; |
---|
51 | % |
---|
52 | % Topography netcdf file name (ETOPO 2) |
---|
53 | % |
---|
54 | topofile='../Topo/etopo2.nc'; |
---|
55 | % |
---|
56 | % GSHSS user defined coastline (see m_map) |
---|
57 | % |
---|
58 | coastfileplot=''; |
---|
59 | % |
---|
60 | % |
---|
61 | %%%%%%%%%%%%%%%%%%% END USERS DEFINED VARIABLES %%%%%%%%%%%%%%%%%%%%%%% |
---|
62 | % |
---|
63 | % Title |
---|
64 | % |
---|
65 | disp(' ') |
---|
66 | disp([' Making the grid: ',grdname]) |
---|
67 | disp(' ') |
---|
68 | disp([' Title: ',title]) |
---|
69 | disp(' ') |
---|
70 | disp([' Resolution: 1/',num2str(1/dl),' deg']) |
---|
71 | % |
---|
72 | % Get the Longitude and latitude |
---|
73 | % |
---|
74 | theta=atan((lat3-lat1)/(lon3-lon1)) |
---|
75 | l=(lon2-lon1)*cos(theta)+(lat2-lat1)*sin(theta) |
---|
76 | L=(lat2-lat1)*cos(theta)-(lon2-lon1)*sin(theta) |
---|
77 | x=(0:dl:l); |
---|
78 | y=(0:dl:L); |
---|
79 | [xr,yr]=meshgrid(x,y); |
---|
80 | Lonr=lon1+xr*cos(theta)-yr*sin(theta); |
---|
81 | Latr=lat1+xr*sin(theta)+yr*cos(theta); |
---|
82 | [Lonu,Lonv,Lonp]=rho2uvp(Lonr); |
---|
83 | [Latu,Latv,Latp]=rho2uvp(Latr); |
---|
84 | % |
---|
85 | % Create the grid file |
---|
86 | % |
---|
87 | disp(' ') |
---|
88 | disp(' Create the grid file...') |
---|
89 | [M,L]=size(Latp); |
---|
90 | disp([' LLm = ',num2str(L-1)]) |
---|
91 | disp([' MMm = ',num2str(M-1)]) |
---|
92 | create_grid(L,M,grdname,title) |
---|
93 | % |
---|
94 | % Fill the grid file |
---|
95 | % |
---|
96 | disp(' ') |
---|
97 | disp(' Fill the grid file...') |
---|
98 | nc=netcdf(grdname,'write'); |
---|
99 | nc{'lat_u'}(:)=Latu; |
---|
100 | nc{'lon_u'}(:)=Lonu; |
---|
101 | nc{'lat_v'}(:)=Latv; |
---|
102 | nc{'lon_v'}(:)=Lonv; |
---|
103 | nc{'lat_rho'}(:)=Latr; |
---|
104 | nc{'lon_rho'}(:)=Lonr; |
---|
105 | nc{'lat_psi'}(:)=Latp; |
---|
106 | nc{'lon_psi'}(:)=Lonp; |
---|
107 | result=close(nc); |
---|
108 | % |
---|
109 | % Compute the metrics |
---|
110 | % |
---|
111 | disp(' ') |
---|
112 | disp(' Compute the metrics...') |
---|
113 | [pm,pn,dndx,dmde]=get_metrics(grdname); |
---|
114 | xr=0.*pm; |
---|
115 | yr=xr; |
---|
116 | for i=1:L |
---|
117 | xr(:,i+1)=xr(:,i)+2./(pm(:,i+1)+pm(:,i)); |
---|
118 | end |
---|
119 | for j=1:M |
---|
120 | yr(j+1,:)=yr(j,:)+2./(pn(j+1,:)+pn(j,:)); |
---|
121 | end |
---|
122 | [xu,xv,xp]=rho2uvp(xr); |
---|
123 | [yu,yv,yp]=rho2uvp(yr); |
---|
124 | dx=1./pm; |
---|
125 | dy=1./pn; |
---|
126 | dxmax=max(max(dx/1000)); |
---|
127 | dxmin=min(min(dx/1000)); |
---|
128 | dymax=max(max(dy/1000)); |
---|
129 | dymin=min(min(dy/1000)); |
---|
130 | disp(' ') |
---|
131 | disp([' Min dx=',num2str(dxmin),' km - Max dx=',num2str(dxmax),' km']) |
---|
132 | disp([' Min dy=',num2str(dymin),' km - Max dy=',num2str(dymax),' km']) |
---|
133 | % |
---|
134 | % Add topography from topofile |
---|
135 | % |
---|
136 | disp(' ') |
---|
137 | disp(' Add topography...') |
---|
138 | h=add_topo(grdname,topofile); |
---|
139 | % |
---|
140 | % Compute the mask |
---|
141 | % |
---|
142 | maskr=h>0; |
---|
143 | maskr=process_mask(maskr); |
---|
144 | if config=='kago' |
---|
145 | maskr(Latr>45.5)=0; |
---|
146 | end |
---|
147 | if config=='test' |
---|
148 | maskr(Lonr<41.5 & Latr<14.4)=0.; |
---|
149 | maskr(Lonr<40.5 & Latr<14.6)=0.; |
---|
150 | end |
---|
151 | [masku,maskv,maskp]=uvp_mask(maskr); |
---|
152 | % |
---|
153 | % Smooth the topography |
---|
154 | % |
---|
155 | h = smoothgrid(h,hmin,rtarget); |
---|
156 | % |
---|
157 | % Angle between XI-axis and the direction |
---|
158 | % to the EAST at RHO-points [radians]. |
---|
159 | % |
---|
160 | angle=get_angle(Latu,Lonu); |
---|
161 | % |
---|
162 | % Coriolis parameter |
---|
163 | % |
---|
164 | f=4*pi*sin(pi*Latr/180)/(24*3600); |
---|
165 | % |
---|
166 | % Write it down |
---|
167 | % |
---|
168 | disp(' ') |
---|
169 | disp(' Write it down...') |
---|
170 | nc=netcdf(grdname,'write'); |
---|
171 | nc{'h'}(:)=h; |
---|
172 | nc{'pm'}(:)=pm; |
---|
173 | nc{'pn'}(:)=pn; |
---|
174 | nc{'dndx'}(:)=dndx; |
---|
175 | nc{'dmde'}(:)=dmde; |
---|
176 | nc{'mask_u'}(:)=masku; |
---|
177 | nc{'mask_v'}(:)=maskv; |
---|
178 | nc{'mask_psi'}(:)=maskp; |
---|
179 | nc{'mask_rho'}(:)=maskr; |
---|
180 | nc{'x_u'}(:)=xu; |
---|
181 | nc{'y_u'}(:)=yu; |
---|
182 | nc{'x_v'}(:)=xv; |
---|
183 | nc{'y_v'}(:)=yv; |
---|
184 | nc{'x_rho'}(:)=xr; |
---|
185 | nc{'y_rho'}(:)=yr; |
---|
186 | nc{'x_psi'}(:)=xp; |
---|
187 | nc{'y_psi'}(:)=yp; |
---|
188 | nc{'angle'}(:)=angle; |
---|
189 | nc{'f'}(:)=f; |
---|
190 | nc{'spherical'}(:)='T'; |
---|
191 | result=close(nc); |
---|
192 | disp(' ') |
---|
193 | disp([' Size of the grid: L = ',... |
---|
194 | num2str(L),' - M = ',num2str(M)]) |
---|
195 | % |
---|
196 | % make a plot |
---|
197 | % |
---|
198 | disp(' ') |
---|
199 | disp(' Do a plot...') |
---|
200 | themask=ones(size(maskr)); |
---|
201 | themask(maskr==0)=NaN; |
---|
202 | domaxis=[min(min(Lonr)) max(max(Lonr)) min(min(Latr)) max(max(Latr))]; |
---|
203 | colaxis=[min(min(h)) max(max(h))]; |
---|
204 | fixcolorbar([0.25 0.05 0.5 0.03],colaxis,... |
---|
205 | 'Topography',10) |
---|
206 | width=1; |
---|
207 | height=0.8; |
---|
208 | subplot('position',[0. 0.14 width height]) |
---|
209 | m_proj('mercator',... |
---|
210 | 'lon',[domaxis(1) domaxis(2)],... |
---|
211 | 'lat',[domaxis(3) domaxis(4)]); |
---|
212 | m_pcolor(Lonr,Latr,h.*themask); |
---|
213 | %shading interp |
---|
214 | caxis(colaxis) |
---|
215 | hold on |
---|
216 | m_contour(Lonr,Latr,h,[hmin 100 200 500 1000 2000 4000],'k--'); |
---|
217 | if ~isempty(coastfileplot) |
---|
218 | m_usercoast(coastfileplot,'patch',[.9 .9 .9]); |
---|
219 | else |
---|
220 | % m_gshhs_l('patch',[.9 .9 .9]) |
---|
221 | end |
---|
222 | m_grid('box','fancy',... |
---|
223 | 'xtick',5,'ytick',5,'tickdir','out',... |
---|
224 | 'fontsize',7); |
---|
225 | hold off |
---|
226 | % |
---|
227 | % End |
---|
228 | % |
---|
229 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|