1 | function geost_currents_bry(bryname,grdname,Zbryname,frcname,zref,obcndx) |
---|
2 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
3 | % |
---|
4 | % compute SSH and the geostrophic currents from Hydrology data |
---|
5 | % |
---|
6 | % |
---|
7 | % Further Information: |
---|
8 | % http://www.brest.ird.fr/Roms_tools/ |
---|
9 | % |
---|
10 | % This file is part of ROMSTOOLS |
---|
11 | % |
---|
12 | % ROMSTOOLS is free software; you can redistribute it and/or modify |
---|
13 | % it under the terms of the GNU General Public License as published |
---|
14 | % by the Free Software Foundation; either version 2 of the License, |
---|
15 | % or (at your option) any later version. |
---|
16 | % |
---|
17 | % ROMSTOOLS is distributed in the hope that it will be useful, but |
---|
18 | % WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
19 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
20 | % GNU General Public License for more details. |
---|
21 | % |
---|
22 | % You should have received a copy of the GNU General Public License |
---|
23 | % along with this program; if not, write to the Free Software |
---|
24 | % Foundation, Inc., 59 Temple Place, Suite 330, Boston, |
---|
25 | % MA 02111-1307 USA |
---|
26 | % |
---|
27 | % Adapted from a previous program of Patrick Marchesiello (IRD). |
---|
28 | % |
---|
29 | % Copyright (c) 2001-2006 by Patrick Marchesiello and Pierrick Penven |
---|
30 | % e-mail:Pierrick.Penven@ird.fr |
---|
31 | % |
---|
32 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
33 | rho0=1025; % Bousinesq background density [kg.m-3] |
---|
34 | g=9.8; % Gravity acceleration [m.s-2] |
---|
35 | De=40; % Ekman depth [m] |
---|
36 | % |
---|
37 | % grid parameters |
---|
38 | % |
---|
39 | % disp(' Read grid parameters ...'); |
---|
40 | nc=netcdf(grdname); |
---|
41 | L=length(nc('xi_rho')); |
---|
42 | M=length(nc('eta_rho')); |
---|
43 | latu=nc{'lat_u'}(:); |
---|
44 | lonu=nc{'lon_u'}(:); |
---|
45 | lonv=nc{'lon_v'}(:); |
---|
46 | latv=nc{'lat_v'}(:); |
---|
47 | lat=nc{'lat_rho'}(:,1); |
---|
48 | lon=nc{'lon_rho'}(1,:); |
---|
49 | |
---|
50 | if obcndx==1 |
---|
51 | h=nc{'h'}(1,:); |
---|
52 | pm=nc{'pm'}(1,:); |
---|
53 | pn=nc{'pn'}(1,:); |
---|
54 | f=nc{'f'}(1,:); |
---|
55 | rmask=nc{'mask_rho'}(1,:); |
---|
56 | umask=nc{'mask_u'}(1,:); |
---|
57 | vmask=nc{'mask_v'}(1,:); |
---|
58 | suffix='_south'; |
---|
59 | elseif obcndx==2 |
---|
60 | h=(nc{'h'}(:,L))'; |
---|
61 | pm=(nc{'pm'}(:,L))'; |
---|
62 | pn=(nc{'pn'}(:,L))'; |
---|
63 | f=(nc{'f'}(:,L))'; |
---|
64 | rmask=(nc{'mask_rho'}(:,L))'; |
---|
65 | umask=(nc{'mask_u'}(:,L-1))'; |
---|
66 | vmask=(nc{'mask_v'}(:,L))'; |
---|
67 | suffix='_east'; |
---|
68 | elseif obcndx==3 |
---|
69 | h=nc{'h'}(M,:); |
---|
70 | pm=nc{'pm'}(M,:); |
---|
71 | pn=nc{'pn'}(M,:); |
---|
72 | f=nc{'f'}(M,:); |
---|
73 | rmask=nc{'mask_rho'}(M,:); |
---|
74 | umask=nc{'mask_u'}(M,:); |
---|
75 | vmask=nc{'mask_v'}(M-1,:); |
---|
76 | suffix='_north'; |
---|
77 | elseif obcndx==4 |
---|
78 | h=(nc{'h'}(:,1))'; |
---|
79 | pm=(nc{'pm'}(:,1))'; |
---|
80 | pn=(nc{'pn'}(:,1))'; |
---|
81 | f=(nc{'f'}(:,1))'; |
---|
82 | rmask=(nc{'mask_rho'}(:,1))'; |
---|
83 | umask=(nc{'mask_u'}(:,1))'; |
---|
84 | vmask=(nc{'mask_v'}(:,1))'; |
---|
85 | suffix='_west'; |
---|
86 | end |
---|
87 | Nx=length(h); |
---|
88 | close(nc) |
---|
89 | % |
---|
90 | % Levitus vertical levels |
---|
91 | % |
---|
92 | noa=netcdf(Zbryname); |
---|
93 | Z=-noa{'Z'}(:); |
---|
94 | NL=length(Z); |
---|
95 | time=noa{'bry_time'}(:); |
---|
96 | tlen=length(time); |
---|
97 | % |
---|
98 | % Model grid vertical levels |
---|
99 | % |
---|
100 | nc=netcdf(bryname,'write'); |
---|
101 | theta_s = nc{'theta_s'}(:); |
---|
102 | theta_b = nc{'theta_b'}(:); |
---|
103 | hc = nc{'hc'}(:); |
---|
104 | N = length(nc('s_rho')); |
---|
105 | % |
---|
106 | % get the reference level |
---|
107 | % |
---|
108 | kref=min(find(Z<=zref)); |
---|
109 | if isempty(kref); |
---|
110 | kref=length(Z); |
---|
111 | disp(['Warning zref not found. Taking :',num2str(Z(kref))]) |
---|
112 | end |
---|
113 | Z=Z(1:kref); |
---|
114 | z=reshape(Z,kref,1,1); |
---|
115 | z=repmat(z,[1 Nx]); |
---|
116 | % |
---|
117 | % Open the forcing file |
---|
118 | % |
---|
119 | if ~isempty(frcname) |
---|
120 | nfrc=netcdf(frcname); |
---|
121 | end |
---|
122 | %%%%%%%%%%%%%%%%%%% |
---|
123 | % START MAIN LOOP |
---|
124 | %%%%%%%%%%%%%%%%%%% |
---|
125 | % |
---|
126 | % loop on time |
---|
127 | % |
---|
128 | for l=1:tlen |
---|
129 | %for l=1:1 |
---|
130 | disp(['time index: ',num2str(l),' of total: ',num2str(tlen)]) |
---|
131 | % |
---|
132 | % read T and S |
---|
133 | % |
---|
134 | T3d=squeeze(noa{['temp',suffix]}(l,1:kref,:)); |
---|
135 | S3d=squeeze(noa{['salt',suffix]}(l,1:kref,:)); |
---|
136 | Ts=squeeze(T3d(1,:)); |
---|
137 | Ss=squeeze(S3d(1,:)); |
---|
138 | rhos=rho_eos(Ts,Ss,0); |
---|
139 | rho=rho_eos(T3d,S3d,z); |
---|
140 | rho_w=.5*(rho(1:kref-1,:)+rho(2:kref,:)); |
---|
141 | z_w=.5*(z(1:kref-1,:)+z(2:kref,:)); |
---|
142 | dz_w=z(1:kref-1,:)-z(2:kref,:); |
---|
143 | % |
---|
144 | % COMPUTE PRESSURE |
---|
145 | % |
---|
146 | % disp('Pressure field and sea level ...') |
---|
147 | pres=0*T3d; |
---|
148 | % initialize pressure at kref in Pascal |
---|
149 | pres(kref,:)=-zref*1.e4; |
---|
150 | for k=kref-1:-1:1; |
---|
151 | pres(k,:)=pres(k+1,:)-rho_w(k,:).*g.*dz_w(k,:); |
---|
152 | end |
---|
153 | % |
---|
154 | % compute SSH |
---|
155 | % |
---|
156 | ssh=(squeeze(pres(1,:))./(rhos*g)); |
---|
157 | % avgssh=sum(rmask.*ssh./(pm.*pn))/sum(rmask./(pm.*pn)); |
---|
158 | % ssh=ssh-avgssh; |
---|
159 | % avgp=squeeze(tridim(avgssh.*rhos*g,kref)); |
---|
160 | % pres=pres-avgp; |
---|
161 | % |
---|
162 | % COMPUTE GEOSTROPHIC BAROCLINIC VELOCITIES |
---|
163 | % |
---|
164 | % disp('Baroclinic geostrophic component ...') |
---|
165 | pn3d=squeeze(tridim(pn,kref)); |
---|
166 | pm3d=squeeze(tridim(pm,kref)); |
---|
167 | f3d=squeeze(tridim(f,kref)); |
---|
168 | m3d=squeeze(tridim(rmask,kref)); |
---|
169 | if obcndx==1 | obcndx==3 |
---|
170 | p_u=0.5*(pres(:,1:Nx-1)+pres(:,2:Nx)); |
---|
171 | px(:,2:Nx-1)=p_u(:,2:Nx-1)-p_u(:,1:Nx-2); |
---|
172 | px(:,1)=2.*px(:,2)-px(:,3); |
---|
173 | px(:,Nx)=2.*px(:,Nx-1)-px(:,Nx-2); |
---|
174 | v_r=m3d.*pm3d.*px./(rho0*f3d); |
---|
175 | u_r=0*v_r; |
---|
176 | end |
---|
177 | if obcndx==2 | obcndx==4 |
---|
178 | p_v=0.5*(pres(:,1:Nx-1)+pres(:,2:Nx)); |
---|
179 | py(:,2:Nx-1)=p_v(:,2:Nx-1)-p_v(:,1:Nx-2); |
---|
180 | py(:,1)=2.*py(:,2)-py(:,3); |
---|
181 | py(:,Nx)=2.*py(:,Nx-1)-py(:,Nx-2); |
---|
182 | u_r=-m3d.*pn3d.*py./(rho0*f3d); |
---|
183 | v_r=0*u_r; |
---|
184 | end |
---|
185 | % |
---|
186 | % Ekman transport |
---|
187 | % |
---|
188 | if ~isempty(frcname) |
---|
189 | % disp('Add the Ekman transport') |
---|
190 | if obcndx==1 |
---|
191 | tmp=squeeze(nfrc{'sustr'}(l,1,:)); |
---|
192 | sustr=0*h; |
---|
193 | sustr(2:end-1)=0.5*(tmp(1:end-1)+tmp(2:end)); |
---|
194 | sustr(1)=sustr(2); |
---|
195 | sustr(end)=sustr(end-1); |
---|
196 | svstr=squeeze(nfrc{'svstr'}(l,1,:)); |
---|
197 | elseif obcndx==2 |
---|
198 | sustr=(squeeze(nfrc{'sustr'}(l,:,L-1)))'; |
---|
199 | svstr=0*h; |
---|
200 | tmp=(squeeze(nfrc{'svstr'}(l,:,L)))'; |
---|
201 | svstr(2:end-1)=0.5*(tmp(1:end-1)+tmp(2:end)); |
---|
202 | svstr(1)=svstr(2); |
---|
203 | svstr(end)=svstr(end-1); |
---|
204 | elseif obcndx==3 |
---|
205 | tmp=squeeze(nfrc{'sustr'}(l,M,:)); |
---|
206 | sustr=0*h; |
---|
207 | sustr(2:end-1)=0.5*(tmp(1:end-1)+tmp(2:end)); |
---|
208 | sustr(1)=sustr(2); |
---|
209 | sustr(end)=sustr(end-1); |
---|
210 | svstr=squeeze(nfrc{'svstr'}(l,M-1,:)); |
---|
211 | elseif obcndx==4 |
---|
212 | sustr=(squeeze(nfrc{'sustr'}(l,:,1)))'; |
---|
213 | svstr=0*h; |
---|
214 | tmp=(squeeze(nfrc{'svstr'}(l,:,1)))'; |
---|
215 | svstr(2:end-1)=0.5*(tmp(1:end-1)+tmp(2:end)); |
---|
216 | svstr(1)=svstr(2); |
---|
217 | svstr(end)=svstr(end-1); |
---|
218 | end |
---|
219 | k_ekman=min(find(Z<=-De)); |
---|
220 | u_r(1:k_ekman,:)=u_r(1:k_ekman,:)+squeeze(tridim(rmask.*... |
---|
221 | svstr./(rho0*De*f),k_ekman)); |
---|
222 | v_r(1:k_ekman,:)=v_r(1:k_ekman,:)+squeeze(tridim(-rmask.*... |
---|
223 | sustr./(rho0*De*f),k_ekman)); |
---|
224 | end |
---|
225 | |
---|
226 | |
---|
227 | |
---|
228 | % Replace/interpolate Equatorial values |
---|
229 | % |
---|
230 | |
---|
231 | if obcndx==2 | obcndx==4 |
---|
232 | |
---|
233 | equatlat=(lat>=-2 & lat<=2); |
---|
234 | if sum(sum(equatlat))~=0 |
---|
235 | % disp('Extrapole values outside the Equator') |
---|
236 | D=find(~equatlat); |
---|
237 | |
---|
238 | |
---|
239 | if length(D)~=0 |
---|
240 | |
---|
241 | for k=1:kref |
---|
242 | u_r(k,:)=interp1(lat(D),u_r(k,D),lat,'spline','extrap'); |
---|
243 | v_r(k,:)=interp1(lat(D),v_r(k,D),lat,'spline','extrap'); |
---|
244 | |
---|
245 | end |
---|
246 | |
---|
247 | else |
---|
248 | disp('No values outside the Equator to extrapole') |
---|
249 | end |
---|
250 | end |
---|
251 | |
---|
252 | elseif obcndx==1 | obcndx==3 |
---|
253 | |
---|
254 | equatlat=(lat>=-2 & lat<=2); |
---|
255 | if sum(sum(equatlat))~=0 |
---|
256 | % disp('Extrapole values outside the Equator') |
---|
257 | D=find(~equatlat); |
---|
258 | |
---|
259 | |
---|
260 | if length(D)~=0 |
---|
261 | |
---|
262 | for k=1:kref |
---|
263 | u_r(k,:)=interp1(lon(D),u_r(k,D),lon,'spline','extrap'); |
---|
264 | v_r(k,:)=interp1(lon(D),v_r(k,D),lon,'spline','extrap'); |
---|
265 | |
---|
266 | end |
---|
267 | |
---|
268 | else |
---|
269 | disp('No values outside the Equator to extrapole') |
---|
270 | end |
---|
271 | end |
---|
272 | |
---|
273 | end |
---|
274 | |
---|
275 | |
---|
276 | % |
---|
277 | % Masking |
---|
278 | % |
---|
279 | umask3d=squeeze(tridim(umask,kref)); |
---|
280 | vmask3d=squeeze(tridim(vmask,kref)); |
---|
281 | if obcndx==1 | obcndx==3 |
---|
282 | u=umask3d.*0.5.*(u_r(:,1:end-1)+u_r(:,2:end)); |
---|
283 | v=vmask3d.*v_r; |
---|
284 | end |
---|
285 | if obcndx==2 | obcndx==4 |
---|
286 | u=umask3d.*u_r; |
---|
287 | v=vmask3d.*0.5.*(v_r(:,1:end-1)+v_r(:,2:end)); |
---|
288 | end |
---|
289 | ssh=ssh.*rmask; |
---|
290 | % |
---|
291 | % Vertical interpolation of baroclinic fields |
---|
292 | % |
---|
293 | % disp('Vertical interpolation ...') |
---|
294 | zroms=squeeze(zlevs(h,0*h,theta_s,theta_b,hc,N,'r')); |
---|
295 | if obcndx==1 | obcndx==3 |
---|
296 | zu=0.5*(zroms(:,1:end-1)+zroms(:,2:end)); |
---|
297 | zv=zroms; |
---|
298 | end |
---|
299 | if obcndx==2 | obcndx==4 |
---|
300 | zu=zroms; |
---|
301 | zv=0.5*(zroms(:,1:end-1)+zroms(:,2:end)); |
---|
302 | end |
---|
303 | % |
---|
304 | % Add non gradient velocities at the top and nul velocities |
---|
305 | % at -10000m for vertical extrapolation. |
---|
306 | % |
---|
307 | u=cat(1,u(1,:),u); |
---|
308 | v=cat(1,v(1,:),v); |
---|
309 | u=flipdim(cat(1,u,0*u(1,:)),1); |
---|
310 | v=flipdim(cat(1,v,0*v(1,:)),1); |
---|
311 | z1=flipud([100;Z;-10000]); |
---|
312 | % |
---|
313 | % Do the interpolation |
---|
314 | % |
---|
315 | u=ztosigma_1d(u,zu,z1); |
---|
316 | v=ztosigma_1d(v,zv,z1); |
---|
317 | % |
---|
318 | % Barotropic velocities |
---|
319 | % |
---|
320 | % disp('Barotropic component ...') |
---|
321 | zw=squeeze(zlevs(h,0*h,theta_s,theta_b,hc,N,'w')); |
---|
322 | dz=zw(2:end,:)-zw(1:end-1,:); |
---|
323 | if obcndx==1 | obcndx==3 |
---|
324 | dzu=0.5*(dz(:,1:end-1)+dz(:,2:end)); |
---|
325 | dzv=dz; |
---|
326 | end |
---|
327 | if obcndx==2 | obcndx==4 |
---|
328 | dzu=dz; |
---|
329 | dzv=0.5*(dz(:,1:end-1)+dz(:,2:end)); |
---|
330 | end |
---|
331 | hu=sum(dzu.*u); |
---|
332 | hv=sum(dzv.*v); |
---|
333 | D_u=sum(dzu); |
---|
334 | D_v=sum(dzv); |
---|
335 | ubar=hu./D_u; |
---|
336 | vbar=hv./D_v; |
---|
337 | % |
---|
338 | % corners (beurk.....) |
---|
339 | % |
---|
340 | ubar(1)=0.5*ubar(2); |
---|
341 | ubar(end)=0.5*ubar(end-1); |
---|
342 | vbar(1)=0.5*vbar(2); |
---|
343 | vbar(end)=0.5*vbar(end-1); |
---|
344 | u(:,1)=0.5*u(:,2); |
---|
345 | u(:,end)=0.5*u(:,end-1); |
---|
346 | v(:,1)=0.5*v(:,2); |
---|
347 | v(:,end)=0.5*v(:,end-1); |
---|
348 | % |
---|
349 | % Write into file |
---|
350 | % |
---|
351 | % disp('Writes into climatology file ...') |
---|
352 | nc{['u',suffix]}(l,:,:)=u; |
---|
353 | nc{['v',suffix]}(l,:,:)=v; |
---|
354 | nc{['ubar',suffix]}(l,:)=ubar; |
---|
355 | nc{['vbar',suffix]}(l,:)=vbar; |
---|
356 | nc{['zeta',suffix]}(l,:)=ssh; |
---|
357 | end |
---|
358 | close(nc) |
---|
359 | % |
---|
360 | % Close the forcing file |
---|
361 | % |
---|
362 | if ~isempty(frcname) |
---|
363 | close(nfrc); |
---|
364 | end |
---|