1 |
function interpickups(dirin,dirout,varargin) |
2 |
% function interpickups(dirin,dirout,snap) |
3 |
% |
4 |
% This function interpolates the data in |
5 |
% a set of mnc pickup files and grid files from the MITgcm |
6 |
% given in dirin/pickup.*.nc and in dirin/grid.*.nc |
7 |
% (this can be a global file or a collection of per-tile pickup files) |
8 |
% to the pickup files dirout/pickup.*.nc and the grid files dirout/grid.*.nc |
9 |
% (this, too, can be a global file or a collection of per-tile pickup files) |
10 |
% |
11 |
% Extrapolation takes place near the domain edges if domain size |
12 |
% of pickout is larger than that of pickin. |
13 |
% |
14 |
% The number of vertical levels must be the same in the two sets of files. |
15 |
% |
16 |
% Snap is an optional argument if there is more than one timestep |
17 |
% in the file. The default is 1. |
18 |
% |
19 |
% May be fishy near boundaries if grid is not uniform... |
20 |
|
21 |
% $Header: $ |
22 |
% $Name: $ |
23 |
|
24 |
if nargin==2 |
25 |
snap=1 |
26 |
else |
27 |
snap=varargin{1} |
28 |
endif |
29 |
|
30 |
if (strcmp(dirin,dirout)) |
31 |
error('dir','You cant use the same input and output directories!') |
32 |
end |
33 |
|
34 |
pickin=dir([dirin '/pickup.*.nc']) |
35 |
gridin=dir([dirin '/grid.*.nc']) |
36 |
if length(pickin)~=length(gridin) |
37 |
error('in','Incompatible number of input pickups and gridfiles') |
38 |
end |
39 |
|
40 |
pickout=dir([dirout '/pickup.*.nc']) |
41 |
gridout=dir([dirout '/grid.*.nc']) |
42 |
if length(pickout)~=length(gridout) |
43 |
error('out','Incompatible number of output pickups and gridfiles') |
44 |
end |
45 |
|
46 |
%%%%%%%%%%%%INPUT SANITY |
47 |
|
48 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
49 |
gin=netcdf([dirin '/' gridin(1).name],'nowrite'); |
50 |
Zcomp=fin{'Z'}(:); |
51 |
gZcomp=gin{'Z'}(:); |
52 |
if (sum(Zcomp~=gZcomp)>0) |
53 |
error('in','Incompatible Z-axis input pickup and gridfile: 1') |
54 |
end |
55 |
|
56 |
Xin=fin{'X'}(:); |
57 |
gXin=gin{'X'}(:); |
58 |
|
59 |
if (sum(gXin~=gXin)>0) |
60 |
error('in','Incompatible x-axis input pickups and gridfile: 1') |
61 |
end |
62 |
|
63 |
Yin=fin{'Y'}(:); |
64 |
gYin=gin{'Y'}(:); |
65 |
if (sum(gYin~=gYin)>0) |
66 |
error('in','Incompatible y-axis input pickups and gridfile: 1') |
67 |
end |
68 |
|
69 |
Xp1in=fin{'Xp1'}(:); |
70 |
gXp1in=gin{'Xp1'}(:); |
71 |
if (sum(gXp1in~=gXp1in)>0) |
72 |
error('in','Incompatible x-axis input pickups and gridfile: 1') |
73 |
end |
74 |
|
75 |
Yp1in=fin{'Yp1'}(:); |
76 |
gYp1in=gin{'Yp1'}(:); |
77 |
if (sum(gYp1in~=gYp1in)>0) |
78 |
error('in','Incompatible y-axis input pickups and gridfile: 1') |
79 |
end |
80 |
|
81 |
close(fin) |
82 |
close(gin) |
83 |
|
84 |
for i=2:length(pickin) |
85 |
fin=netcdf([dirin '/' pickin(i).name],'nowrite'); |
86 |
Z=fin{'Z'}(:); |
87 |
if (sum(Zcomp~=Z)>0) |
88 |
error('Z','Incompatible vertical axes in input pickups:',num2str(i)) |
89 |
end |
90 |
|
91 |
gin=netcdf([dirin '/' pickin(i).name],'nowrite'); |
92 |
Z=gin{'Z'}(:); |
93 |
if (sum(Zcomp~=Z)>0) |
94 |
error('Z','Incompatible vertical axes in input gridfiles:',num2str(i)) |
95 |
end |
96 |
|
97 |
Xin=sort([Xin;fin{'X'}(:)]); |
98 |
Xp1in=sort([Xp1in;fin{'Xp1'}(:)]); |
99 |
|
100 |
gXin=sort([gXin;gin{'X'}(:)]); |
101 |
gXp1in=sort([gXp1in;gin{'Xp1'}(:)]); |
102 |
|
103 |
if (sum(gXin~=Xin)>0) |
104 |
error('X','Incompatible x-axes in input files:',num2str(i)) |
105 |
end |
106 |
|
107 |
Yin=sort([Yin;fin{'Y'}(:)]); |
108 |
Yp1in=sort([Yp1in;fin{'Yp1'}(:)]); |
109 |
|
110 |
gYin=sort([gYin;fin{'Y'}(:)]); |
111 |
gYp1in=sort([gYp1in;fin{'Yp1'}(:)]); |
112 |
|
113 |
if (sum(gYin~=Yin)>0) |
114 |
error('Y','Incompatible y-axes in input files:',num2str(i)) |
115 |
end |
116 |
|
117 |
close(fin); |
118 |
close(gin); |
119 |
end |
120 |
|
121 |
store=[Xin(1)]; |
122 |
for i=2:length(Xin) |
123 |
if Xin(i-1)~=Xin(i) |
124 |
store(end+1)=Xin(i); |
125 |
end |
126 |
end |
127 |
Xin=store'; |
128 |
clear gXin |
129 |
|
130 |
store=[Xp1in(1)]; |
131 |
for i=2:length(Xp1in) |
132 |
if Xp1in(i-1)~=Xp1in(i) |
133 |
store(end+1)=Xp1in(i); |
134 |
end |
135 |
end |
136 |
Xp1in=store'; |
137 |
clear gXp1in |
138 |
|
139 |
store=[Yin(1)]; |
140 |
for i=2:length(Yin) |
141 |
if Yin(i-1)~=Yin(i) |
142 |
store(end+1)=Yin(i); |
143 |
end |
144 |
end |
145 |
Yin=store'; |
146 |
clear gYin |
147 |
|
148 |
store=[Yp1in(1)]; |
149 |
for i=2:length(Yp1in) |
150 |
if Yp1in(i-1)~=Yp1in(i) |
151 |
store(end+1)=Yp1in(i); |
152 |
end |
153 |
end |
154 |
Yp1in=store'; |
155 |
clear gYp1in |
156 |
|
157 |
%%%%%%%%%%%%%%% OUTPUT SANITY |
158 |
fout=netcdf([dirout '/' pickout(1).name],'nowrite'); |
159 |
gout=netcdf([dirout '/' gridout(1).name],'nowrite'); |
160 |
Zcomp=fout{'Z'}(:); |
161 |
gZcomp=gout{'Z'}(:); |
162 |
if (sum(Zcomp~=gZcomp)>0) |
163 |
error('out','Incompatible Z-axis output pickup and gridfile: 1') |
164 |
end |
165 |
|
166 |
Xout=fout{'X'}(:); |
167 |
gXout=gout{'X'}(:); |
168 |
|
169 |
if (sum(gXout~=gXout)>0) |
170 |
error('out','Incompatible x-axis output pickups and gridfile: 1') |
171 |
end |
172 |
|
173 |
Yout=fout{'Y'}(:); |
174 |
gYout=gout{'Y'}(:); |
175 |
if (sum(gYout~=gYout)>0) |
176 |
error('out','Incompatible y-axis output pickups and gridfile: 1') |
177 |
end |
178 |
|
179 |
Xp1out=fout{'Xp1'}(:); |
180 |
gXp1out=gout{'Xp1'}(:); |
181 |
if (sum(gXp1out~=gXp1out)>0) |
182 |
error('out','Incompatible x-axis output pickups and gridfile: 1') |
183 |
end |
184 |
|
185 |
Yp1out=fout{'Yp1'}(:); |
186 |
gYp1out=gout{'Yp1'}(:); |
187 |
if (sum(gYp1out~=gYp1out)>0) |
188 |
error('out','Incompatible y-axis output pickups and gridfile: 1') |
189 |
end |
190 |
|
191 |
close(fout) |
192 |
close(gout) |
193 |
|
194 |
for i=2:length(pickout) |
195 |
fout=netcdf([dirout '/' pickout(i).name],'nowrite'); |
196 |
Z=fout{'Z'}(:); |
197 |
if (sum(Zcomp~=Z)>0) |
198 |
error('Z','Incompatible vertical axes in output pickups:',num2str(i)) |
199 |
end |
200 |
|
201 |
gout=netcdf([dirout '/' pickout(i).name],'nowrite'); |
202 |
Z=gout{'Z'}(:); |
203 |
if (sum(Zcomp~=Z)>0) |
204 |
error('Z','Incompatible vertical axes in output gridfiles:',num2str(i)) |
205 |
end |
206 |
|
207 |
Xout=sort([Xout;fout{'X'}(:)]); |
208 |
Xp1out=sort([Xp1out;fout{'Xp1'}(:)]); |
209 |
|
210 |
gXout=sort([gXout;gout{'X'}(:)]); |
211 |
gXp1out=sort([gXp1out;gout{'Xp1'}(:)]); |
212 |
|
213 |
if (sum(gXout~=Xout)>0) |
214 |
error('X','Incompatible x-axes in output files:',num2str(i)) |
215 |
end |
216 |
|
217 |
Yout=sort([Yout;fout{'Y'}(:)]); |
218 |
Yp1out=sort([Yp1out;fout{'Yp1'}(:)]); |
219 |
|
220 |
gYout=sort([gYout;fout{'Y'}(:)]); |
221 |
gYp1out=sort([gYp1out;fout{'Yp1'}(:)]); |
222 |
|
223 |
if (sum(gYout~=Yout)>0) |
224 |
error('Y','Incompatible y-axes in output files:',num2str(i)) |
225 |
end |
226 |
|
227 |
close(fout); |
228 |
close(gout); |
229 |
end |
230 |
|
231 |
store=[Xout(1)]; |
232 |
for i=2:length(Xout) |
233 |
if Xout(i-1)~=Xout(i) |
234 |
store(end+1)=Xout(i); |
235 |
end |
236 |
end |
237 |
Xout=store'; |
238 |
clear gXout |
239 |
|
240 |
store=[Xp1out(1)]; |
241 |
for i=2:length(Xp1out) |
242 |
if Xp1out(i-1)~=Xp1out(i) |
243 |
store(end+1)=Xp1out(i); |
244 |
end |
245 |
end |
246 |
Xp1out=store'; |
247 |
clear gXp1out |
248 |
|
249 |
store=[Yout(1)]; |
250 |
for i=2:length(Yout) |
251 |
if Yout(i-1)~=Yout(i) |
252 |
store(end+1)=Yout(i); |
253 |
end |
254 |
end |
255 |
Yout=store'; |
256 |
clear gYout |
257 |
|
258 |
store=[Yp1out(1)]; |
259 |
for i=2:length(Yp1out) |
260 |
if Yp1out(i-1)~=Yp1out(i) |
261 |
store(end+1)=Yp1out(i); |
262 |
end |
263 |
end |
264 |
Yp1out=store'; |
265 |
clear gYp1out |
266 |
|
267 |
|
268 |
% First, do the Centered variables |
269 |
|
270 |
% We assume periodicity as MITgcm does |
271 |
Xinext=[2*Xin(1)-Xin(2);Xin;2*Xin(end)-Xin(end-1)]; |
272 |
Yinext=[2*Yin(1)-Yin(2);Yin;2*Yin(end)-Yin(end-1)]; |
273 |
|
274 |
[xcin,ycin]=meshgrid(Xinext,Yinext); |
275 |
[xcout,ycout]=meshgrid(Xout,Yout); |
276 |
|
277 |
%%%%%%%%%%%%% HFacCoutk |
278 |
|
279 |
HFacoutk=ones([length(Zcomp) size(xcout)]); |
280 |
|
281 |
disp(['Calculating HFacC...']) |
282 |
for j=1:length(gridout) |
283 |
gout=netcdf([dirout '/' gridout(j).name],'nowrite'); |
284 |
Xhere=gout{'X'}(:); |
285 |
Yhere=gout{'Y'}(:); |
286 |
HFacouthere=gout{'HFacC'}(:,:,:); |
287 |
for ii=1:length(Xhere) |
288 |
for jj=1:length(Yhere) |
289 |
[jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout)); |
290 |
HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii); |
291 |
end |
292 |
end |
293 |
close(gout) |
294 |
end |
295 |
clear HFacouthere |
296 |
disp(['Done.']) |
297 |
|
298 |
%%%%%%%%%%%%% ETA |
299 |
|
300 |
Fieldin=NaN*ones(size(xcin)); |
301 |
Fieldout=NaN*ones(size(xcout)); |
302 |
|
303 |
Fieldin=NaN*ones(size(xcin)); |
304 |
Fieldout=NaN*ones(size(xcout)); |
305 |
for j=1:length(pickin) |
306 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
307 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
308 |
Xhere=fin{'X'}(:); |
309 |
Yhere=fin{'Y'}(:); |
310 |
Fieldinhere=fin{'Eta'}(:,:); |
311 |
HFacinhere=squeeze(gin{'HFacC'}(1,:,:)); |
312 |
for ii=1:length(Xhere) |
313 |
for jj=1:length(Yhere) |
314 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
315 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
316 |
if (HFacinhere(jj,ii)==0) |
317 |
Fieldin(jjj,iii)=NaN; |
318 |
end |
319 |
end |
320 |
end |
321 |
close(fin) |
322 |
close(gin) |
323 |
end |
324 |
% This utilizes periodic geometry |
325 |
Fieldin(:,1)=Fieldin(:,end-1); |
326 |
Fieldin(:,end)=Fieldin(:,2); |
327 |
|
328 |
Fieldin(1,:)=Fieldin(end-1,:); |
329 |
Fieldin(end,:)=Fieldin(2,:); |
330 |
|
331 |
disp(['Interpolating Eta ...']) |
332 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
333 |
Fieldout=inpaint_nans(Field0,0); |
334 |
disp('Done.') |
335 |
|
336 |
disp(['Zeroing Eta within topography']) |
337 |
Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:)); |
338 |
disp('Done.') |
339 |
|
340 |
for j=1:length(pickout) |
341 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
342 |
Xhere=fout{'X'}(:); |
343 |
Yhere=fout{'Y'}(:); |
344 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
345 |
fout{'Eta'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
346 |
close(fout) |
347 |
end |
348 |
|
349 |
%%%%%%%%%%%%% EtaH |
350 |
|
351 |
Fieldin=NaN*ones(size(xcin)); |
352 |
Fieldout=NaN*ones(size(xcout)); |
353 |
|
354 |
Fieldin=NaN*ones(size(xcin)); |
355 |
Fieldout=NaN*ones(size(xcout)); |
356 |
for j=1:length(pickin) |
357 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
358 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
359 |
Xhere=fin{'X'}(:); |
360 |
Yhere=fin{'Y'}(:); |
361 |
Fieldinhere=fin{'EtaH'}(:,:); |
362 |
HFacinhere=squeeze(gin{'HFacC'}(1,:,:)); |
363 |
for ii=1:length(Xhere) |
364 |
for jj=1:length(Yhere) |
365 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
366 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
367 |
if (HFacinhere(jj,ii)==0) |
368 |
Fieldin(jjj,iii)=NaN; |
369 |
end |
370 |
end |
371 |
end |
372 |
close(fin) |
373 |
close(gin) |
374 |
end |
375 |
% This utilizes periodic geometry |
376 |
Fieldin(:,1)=Fieldin(:,end-1); |
377 |
Fieldin(:,end)=Fieldin(:,2); |
378 |
|
379 |
Fieldin(1,:)=Fieldin(end-1,:); |
380 |
Fieldin(end,:)=Fieldin(2,:); |
381 |
|
382 |
disp(['Interpolating EtaH ...']) |
383 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
384 |
Fieldout=inpaint_nans(Field0,0); |
385 |
disp('Done.') |
386 |
|
387 |
disp(['Zeroing EtaH within topography']) |
388 |
Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:)); |
389 |
disp('Done.') |
390 |
|
391 |
for j=1:length(pickout) |
392 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
393 |
Xhere=fout{'X'}(:); |
394 |
Yhere=fout{'Y'}(:); |
395 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
396 |
fout{'EtaH'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
397 |
close(fout) |
398 |
end |
399 |
|
400 |
%%%%%%%%%%%%% EtaHdt |
401 |
|
402 |
Fieldin=NaN*ones(size(xcin)); |
403 |
Fieldout=NaN*ones(size(xcout)); |
404 |
|
405 |
Fieldin=NaN*ones(size(xcin)); |
406 |
Fieldout=NaN*ones(size(xcout)); |
407 |
for j=1:length(pickin) |
408 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
409 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
410 |
Xhere=fin{'X'}(:); |
411 |
Yhere=fin{'Y'}(:); |
412 |
Fieldinhere=fin{'dEtaHdt'}(:,:); |
413 |
HFacinhere=squeeze(gin{'HFacC'}(1,:,:)); |
414 |
for ii=1:length(Xhere) |
415 |
for jj=1:length(Yhere) |
416 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
417 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
418 |
if (HFacinhere(jj,ii)==0) |
419 |
Fieldin(jjj,iii)=NaN; |
420 |
end |
421 |
end |
422 |
end |
423 |
close(fin) |
424 |
close(gin) |
425 |
end |
426 |
% This utilizes periodic geometry |
427 |
Fieldin(:,1)=Fieldin(:,end-1); |
428 |
Fieldin(:,end)=Fieldin(:,2); |
429 |
|
430 |
Fieldin(1,:)=Fieldin(end-1,:); |
431 |
Fieldin(end,:)=Fieldin(2,:); |
432 |
|
433 |
disp(['Interpolating dEtaHdt ...']) |
434 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
435 |
Fieldout=inpaint_nans(Field0,0); |
436 |
disp('Done.') |
437 |
|
438 |
disp(['Zeroing dEtaHdt within topography']) |
439 |
Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:)); |
440 |
disp('Done.') |
441 |
|
442 |
for j=1:length(pickout) |
443 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
444 |
Xhere=fout{'X'}(:); |
445 |
Yhere=fout{'Y'}(:); |
446 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
447 |
fout{'dEtaHdt'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
448 |
close(fout) |
449 |
end |
450 |
|
451 |
%S,gSnm1,Temp,gTnm1,phi_nh are on Xc,Rc |
452 |
|
453 |
%%%%%%%%%%%%%% S |
454 |
|
455 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
456 |
for k=1:length(Zcomp) |
457 |
clear Fieldinhere |
458 |
Fieldin=NaN*ones(size(xcin)); |
459 |
Fieldout=NaN*ones(size(xcout)); |
460 |
|
461 |
for j=1:length(pickin) |
462 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
463 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
464 |
Xhere=fin{'X'}(:); |
465 |
Yhere=fin{'Y'}(:); |
466 |
Fieldinhere=squeeze(fin{'S'}(snap,k,:,:)); |
467 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
468 |
for ii=1:length(Xhere) |
469 |
for jj=1:length(Yhere) |
470 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
471 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
472 |
if (HFacinhere(jj,ii)==0) |
473 |
Fieldin(jjj,iii)=NaN; |
474 |
end |
475 |
end |
476 |
end |
477 |
close(fin) |
478 |
close(gin) |
479 |
end |
480 |
% This utilizes periodic geometry |
481 |
Fieldin(:,1)=Fieldin(:,end-1); |
482 |
Fieldin(:,end)=Fieldin(:,2); |
483 |
|
484 |
Fieldin(1,:)=Fieldin(end-1,:); |
485 |
Fieldin(end,:)=Fieldin(2,:); |
486 |
|
487 |
disp(['Interpolating S:',num2str(k),'...']) |
488 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
489 |
Fieldout=inpaint_nans(Field0,0); |
490 |
disp('Done.') |
491 |
|
492 |
disp(['Zeroing S:',num2str(k),' within topography']) |
493 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
494 |
disp('Done.') |
495 |
end |
496 |
|
497 |
for j=1:length(pickout) |
498 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
499 |
Xhere=fout{'X'}(:); |
500 |
Yhere=fout{'Y'}(:); |
501 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
502 |
fout{'S'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
503 |
close(fout) |
504 |
end |
505 |
|
506 |
%%%%%%%%%%%%%% gSnm1 |
507 |
|
508 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
509 |
for k=1:length(Zcomp) |
510 |
Fieldin=NaN*ones(size(xcin)); |
511 |
Fieldout=NaN*ones(size(xcout)); |
512 |
|
513 |
for j=1:length(pickin) |
514 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
515 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
516 |
Xhere=fin{'X'}(:); |
517 |
Yhere=fin{'Y'}(:); |
518 |
Fieldinhere=squeeze(fin{'gSnm1'}(snap,k,:,:)); |
519 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
520 |
for ii=1:length(Xhere) |
521 |
for jj=1:length(Yhere) |
522 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
523 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
524 |
if (HFacinhere(jj,ii)==0) |
525 |
Fieldin(jjj,iii)=NaN; |
526 |
end |
527 |
end |
528 |
end |
529 |
close(fin) |
530 |
close(gin) |
531 |
end |
532 |
% This utilizes periodic geometry |
533 |
Fieldin(:,1)=Fieldin(:,end-1); |
534 |
Fieldin(:,end)=Fieldin(:,2); |
535 |
|
536 |
Fieldin(1,:)=Fieldin(end-1,:); |
537 |
Fieldin(end,:)=Fieldin(2,:); |
538 |
|
539 |
disp(['Interpolating gSnm1:',num2str(k),'...']) |
540 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
541 |
Fieldout=inpaint_nans(Field0,0); |
542 |
disp('Done.') |
543 |
|
544 |
disp(['Zeroing gSnm1:',num2str(k),' within topography']) |
545 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
546 |
disp('Done.') |
547 |
end |
548 |
|
549 |
for j=1:length(pickout) |
550 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
551 |
Xhere=fout{'X'}(:); |
552 |
Yhere=fout{'Y'}(:); |
553 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
554 |
fout{'gSnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
555 |
close(fout) |
556 |
end |
557 |
|
558 |
%%%%%%%%%%%%%% Temp |
559 |
|
560 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
561 |
for k=1:length(Zcomp) |
562 |
Fieldin=NaN*ones(size(xcin)); |
563 |
Fieldout=NaN*ones(size(xcout)); |
564 |
|
565 |
for j=1:length(pickin) |
566 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
567 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
568 |
Xhere=fin{'X'}(:); |
569 |
Yhere=fin{'Y'}(:); |
570 |
Fieldinhere=squeeze(fin{'Temp'}(snap,k,:,:)); |
571 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
572 |
for ii=1:length(Xhere) |
573 |
for jj=1:length(Yhere) |
574 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
575 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
576 |
if (HFacinhere(jj,ii)==0) |
577 |
Fieldin(jjj,iii)=NaN; |
578 |
end |
579 |
end |
580 |
end |
581 |
close(fin) |
582 |
close(gin) |
583 |
end |
584 |
% This utilizes periodic geometry |
585 |
Fieldin(:,1)=Fieldin(:,end-1); |
586 |
Fieldin(:,end)=Fieldin(:,2); |
587 |
|
588 |
Fieldin(1,:)=Fieldin(end-1,:); |
589 |
Fieldin(end,:)=Fieldin(2,:); |
590 |
|
591 |
disp(['Interpolating Temp:',num2str(k),'...']) |
592 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
593 |
Fieldout=inpaint_nans(Field0,0); |
594 |
disp('Done.') |
595 |
|
596 |
disp(['Zeroing Temp:',num2str(k),' within topography']) |
597 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
598 |
disp('Done.') |
599 |
end |
600 |
|
601 |
for j=1:length(pickout) |
602 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
603 |
Xhere=fout{'X'}(:); |
604 |
Yhere=fout{'Y'}(:); |
605 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
606 |
fout{'Temp'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
607 |
close(fout) |
608 |
end |
609 |
|
610 |
%%%%%%%%%%%%%% gTnm1 |
611 |
|
612 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
613 |
for k=1:length(Zcomp) |
614 |
Fieldin=NaN*ones(size(xcin)); |
615 |
Fieldout=NaN*ones(size(xcout)); |
616 |
|
617 |
for j=1:length(pickin) |
618 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
619 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
620 |
Xhere=fin{'X'}(:); |
621 |
Yhere=fin{'Y'}(:); |
622 |
Fieldinhere=squeeze(fin{'gTnm1'}(snap,k,:,:)); |
623 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
624 |
for ii=1:length(Xhere) |
625 |
for jj=1:length(Yhere) |
626 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
627 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
628 |
if (HFacinhere(jj,ii)==0) |
629 |
Fieldin(jjj,iii)=NaN; |
630 |
end |
631 |
end |
632 |
end |
633 |
close(fin) |
634 |
close(gin) |
635 |
end |
636 |
% This utilizes periodic geometry |
637 |
Fieldin(:,1)=Fieldin(:,end-1); |
638 |
Fieldin(:,end)=Fieldin(:,2); |
639 |
|
640 |
Fieldin(1,:)=Fieldin(end-1,:); |
641 |
Fieldin(end,:)=Fieldin(2,:); |
642 |
|
643 |
disp(['Interpolating gTnm1:',num2str(k),'...']) |
644 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
645 |
Fieldout=inpaint_nans(Field0,0); |
646 |
disp('Done.') |
647 |
|
648 |
disp(['Zeroing gTnm1:',num2str(k),' within topography']) |
649 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
650 |
disp('Done.') |
651 |
end |
652 |
|
653 |
for j=1:length(pickout) |
654 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
655 |
Xhere=fout{'X'}(:); |
656 |
Yhere=fout{'Y'}(:); |
657 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
658 |
fout{'gTnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
659 |
close(fout) |
660 |
end |
661 |
|
662 |
%%%%%%%%%%%%%% phi_nh |
663 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
664 |
status=fin{'phi_nh'}; |
665 |
if ~isempty(status) |
666 |
|
667 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
668 |
for k=1:length(Zcomp) |
669 |
Fieldin=NaN*ones(size(xcin)); |
670 |
Fieldout=NaN*ones(size(xcout)); |
671 |
|
672 |
for j=1:length(pickin) |
673 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
674 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
675 |
Xhere=fin{'X'}(:); |
676 |
Yhere=fin{'Y'}(:); |
677 |
Fieldinhere=squeeze(fin{'phi_nh'}(snap,k,:,:)); |
678 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
679 |
for ii=1:length(Xhere) |
680 |
for jj=1:length(Yhere) |
681 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
682 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
683 |
if (HFacinhere(jj,ii)==0) |
684 |
Fieldin(jjj,iii)=NaN; |
685 |
end |
686 |
end |
687 |
end |
688 |
close(fin) |
689 |
close(gin) |
690 |
end |
691 |
% This utilizes periodic geometry |
692 |
Fieldin(:,1)=Fieldin(:,end-1); |
693 |
Fieldin(:,end)=Fieldin(:,2); |
694 |
|
695 |
Fieldin(1,:)=Fieldin(end-1,:); |
696 |
Fieldin(end,:)=Fieldin(2,:); |
697 |
|
698 |
disp(['Interpolating phi_nh:',num2str(k),'...']) |
699 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
700 |
Fieldout=inpaint_nans(Field0,0); |
701 |
disp('Done.') |
702 |
|
703 |
disp(['Zeroing phi_nh:',num2str(k),' within topography']) |
704 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
705 |
disp('Done.') |
706 |
end |
707 |
|
708 |
for j=1:length(pickout) |
709 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
710 |
Xhere=fout{'X'}(:); |
711 |
Yhere=fout{'Y'}(:); |
712 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
713 |
fout{'phi_nh'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
714 |
close(fout) |
715 |
end |
716 |
|
717 |
end %isempty(status) |
718 |
|
719 |
%%%%%%%%%%%%%%%%%% gW |
720 |
%BFK Not sure this is right, since HFacC isn't on the right grid for gW |
721 |
|
722 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
723 |
status=fin{'gW'}; |
724 |
|
725 |
if ~isempty(status) |
726 |
|
727 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
728 |
for k=1:length(Zcomp) |
729 |
Fieldin=NaN*ones(size(xcin)); |
730 |
Fieldout=NaN*ones(size(xcout)); |
731 |
|
732 |
for j=1:length(pickin) |
733 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
734 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
735 |
Xhere=fin{'X'}(:); |
736 |
Yhere=fin{'Y'}(:); |
737 |
Fieldinhere=squeeze(fin{'gW'}(snap,k,:,:)); |
738 |
HFacinhere=squeeze(gin{'HFacC'}(k,:,:)); |
739 |
for ii=1:length(Xhere) |
740 |
for jj=1:length(Yhere) |
741 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
742 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
743 |
if (HFacinhere(jj,ii)==0) |
744 |
Fieldin(jjj,iii)=NaN; |
745 |
end |
746 |
end |
747 |
end |
748 |
close(fin) |
749 |
close(gin) |
750 |
end |
751 |
% This utilizes periodic geometry |
752 |
Fieldin(:,1)=Fieldin(:,end-1); |
753 |
Fieldin(:,end)=Fieldin(:,2); |
754 |
|
755 |
Fieldin(1,:)=Fieldin(end-1,:); |
756 |
Fieldin(end,:)=Fieldin(2,:); |
757 |
|
758 |
disp(['Interpolating gW:',num2str(k),'...']) |
759 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
760 |
Fieldout=inpaint_nans(Field0,0); |
761 |
disp('Done.') |
762 |
|
763 |
disp(['Zeroing gW:',num2str(k),' within topography']) |
764 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
765 |
disp('Done.') |
766 |
end |
767 |
|
768 |
for j=1:length(pickout) |
769 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
770 |
Xhere=fout{'X'}(:); |
771 |
Yhere=fout{'Y'}(:); |
772 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
773 |
fout{'gW'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
774 |
close(fout) |
775 |
end |
776 |
end %isempty(status) |
777 |
|
778 |
%U,gUnm1 is on XU,Rc |
779 |
%Note that there is already a periodic repeat in Uin... |
780 |
Xp1inext=[2*Xp1in(1)-Xp1in(2);Xp1in;2*Xp1in(end)-Xp1in(end-1)]; |
781 |
|
782 |
[xcin,ycin]=meshgrid(Xp1inext,Yinext); |
783 |
[xcout,ycout]=meshgrid(Xp1out,Yout); |
784 |
|
785 |
%%%%%%%%%%%%% HFacWoutk |
786 |
|
787 |
HFacout=ones(size(xcout)); |
788 |
HFacoutk=ones([length(Zcomp) size(xcout)]); |
789 |
|
790 |
disp(['Calculating HFacW...']) |
791 |
for j=1:length(gridout) |
792 |
gout=netcdf([dirout '/' gridout(j).name],'nowrite'); |
793 |
Xhere=gout{'Xp1'}(:); |
794 |
Yhere=gout{'Y'}(:); |
795 |
HFacouthere=gout{'HFacW'}(:,:,:); |
796 |
for ii=1:length(Xhere) |
797 |
for jj=1:length(Yhere) |
798 |
[jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout)); |
799 |
HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii); |
800 |
end |
801 |
end |
802 |
close(gout) |
803 |
end |
804 |
clear HFacouthere |
805 |
disp(['Done.']) |
806 |
|
807 |
%%%%%%%%%%%%%% U |
808 |
|
809 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
810 |
for k=1:length(Zcomp) |
811 |
Fieldin=NaN*ones(size(xcin)); |
812 |
Fieldout=NaN*ones(size(xcout)); |
813 |
|
814 |
for j=1:length(pickin) |
815 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
816 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
817 |
Xhere=fin{'Xp1'}(:); |
818 |
Yhere=fin{'Y'}(:); |
819 |
Fieldinhere=squeeze(fin{'U'}(snap,k,:,:)); |
820 |
HFacinhere=squeeze(gin{'HFacW'}(k,:,:)); |
821 |
for ii=1:length(Xhere) |
822 |
for jj=1:length(Yhere) |
823 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
824 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
825 |
if (HFacinhere(jj,ii)==0) |
826 |
Fieldin(jjj,iii)=NaN; |
827 |
end |
828 |
end |
829 |
end |
830 |
close(fin) |
831 |
close(gin) |
832 |
end |
833 |
% This utilizes periodic geometry |
834 |
%Note that there is already a periodic repeat in U... |
835 |
Fieldin(:,1)=Fieldin(:,end-2); |
836 |
Fieldin(:,end)=Fieldin(:,3); |
837 |
|
838 |
Fieldin(1,:)=Fieldin(end-1,:); |
839 |
Fieldin(end,:)=Fieldin(2,:); |
840 |
|
841 |
disp(['Interpolating U:',num2str(k),'...']) |
842 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
843 |
Fieldout=inpaint_nans(Field0,0); |
844 |
disp('Done.') |
845 |
|
846 |
disp(['Zeroing U:',num2str(k),' within topography']) |
847 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
848 |
disp('Done.') |
849 |
end |
850 |
|
851 |
for j=1:length(pickout) |
852 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
853 |
Xhere=fout{'Xp1'}(:); |
854 |
Yhere=fout{'Y'}(:); |
855 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
856 |
fout{'U'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
857 |
close(fout) |
858 |
end |
859 |
|
860 |
%%%%%%%%%%%%%% gUnm1 |
861 |
|
862 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
863 |
for k=1:length(Zcomp) |
864 |
Fieldin=NaN*ones(size(xcin)); |
865 |
Fieldout=NaN*ones(size(xcout)); |
866 |
|
867 |
for j=1:length(pickin) |
868 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
869 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
870 |
Xhere=fin{'Xp1'}(:); |
871 |
Yhere=fin{'Y'}(:); |
872 |
Fieldinhere=squeeze(fin{'gUnm1'}(snap,k,:,:)); |
873 |
HFacinhere=squeeze(gin{'HFacW'}(k,:,:)); |
874 |
for ii=1:length(Xhere) |
875 |
for jj=1:length(Yhere) |
876 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
877 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
878 |
if (HFacinhere(jj,ii)==0) |
879 |
Fieldin(jjj,iii)=NaN; |
880 |
end |
881 |
end |
882 |
end |
883 |
close(fin) |
884 |
close(gin) |
885 |
end |
886 |
% This utilizes periodic geometry |
887 |
%Note that there is already a periodic repeat in gUnm1... |
888 |
Fieldin(:,1)=Fieldin(:,end-2); |
889 |
Fieldin(:,end)=Fieldin(:,3); |
890 |
|
891 |
Fieldin(1,:)=Fieldin(end-1,:); |
892 |
Fieldin(end,:)=Fieldin(2,:); |
893 |
|
894 |
disp(['Interpolating gUnm1:',num2str(k),'...']) |
895 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
896 |
Fieldout=inpaint_nans(Field0,0); |
897 |
disp('Done.') |
898 |
|
899 |
disp(['Zeroing gUnm1:',num2str(k),' within topography']) |
900 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
901 |
disp('Done.') |
902 |
end |
903 |
|
904 |
for j=1:length(pickout) |
905 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
906 |
Xhere=fout{'Xp1'}(:); |
907 |
Yhere=fout{'Y'}(:); |
908 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
909 |
fout{'gUnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
910 |
close(fout) |
911 |
end |
912 |
|
913 |
%V,gVnm1 is on XV,Rc |
914 |
|
915 |
%Note that there is already a periodic repeat in Vin... |
916 |
Yp1inext=[2*Yp1in(1)-Yp1in(2);Yp1in;2*Yp1in(end)-Yp1in(end-1)]; |
917 |
|
918 |
[xcin,ycin]=meshgrid(Xinext,Yp1inext); |
919 |
[xcout,ycout]=meshgrid(Xout,Yp1out); |
920 |
|
921 |
%%%%%%%%%%%%% HFacSoutk |
922 |
|
923 |
HFacout=ones(size(xcout)); |
924 |
HFacoutk=ones([length(Zcomp) size(xcout)]); |
925 |
|
926 |
disp(['Calculating HFacS...']) |
927 |
for j=1:length(gridout) |
928 |
gout=netcdf([dirout '/' gridout(j).name],'nowrite'); |
929 |
Xhere=gout{'X'}(:); |
930 |
Yhere=gout{'Yp1'}(:); |
931 |
HFacouthere=gout{'HFacS'}(:,:,:); |
932 |
for ii=1:length(Xhere) |
933 |
for jj=1:length(Yhere) |
934 |
[jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout)); |
935 |
HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii); |
936 |
end |
937 |
end |
938 |
close(gout) |
939 |
end |
940 |
clear HFacouthere |
941 |
disp(['Done.']) |
942 |
|
943 |
%%%%%%%%%%%%%% V |
944 |
|
945 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
946 |
for k=1:length(Zcomp) |
947 |
Fieldin=NaN*ones(size(xcin)); |
948 |
Fieldout=NaN*ones(size(xcout)); |
949 |
|
950 |
for j=1:length(pickin) |
951 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
952 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
953 |
Xhere=fin{'X'}(:); |
954 |
Yhere=fin{'Yp1'}(:); |
955 |
Fieldinhere=squeeze(fin{'V'}(snap,k,:,:)); |
956 |
HFacinhere=squeeze(gin{'HFacS'}(k,:,:)); |
957 |
for ii=1:length(Xhere) |
958 |
for jj=1:length(Yhere) |
959 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
960 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
961 |
if (HFacinhere(jj,ii)==0) |
962 |
Fieldin(jjj,iii)=NaN; |
963 |
end |
964 |
end |
965 |
end |
966 |
close(fin) |
967 |
close(gin) |
968 |
end |
969 |
% This utilizes periodic geometry |
970 |
%Note that there is already a periodic repeat in V... |
971 |
Fieldin(:,1)=Fieldin(:,end-1); |
972 |
Fieldin(:,end)=Fieldin(:,2); |
973 |
|
974 |
Fieldin(1,:)=Fieldin(end-2,:); |
975 |
Fieldin(end,:)=Fieldin(3,:); |
976 |
|
977 |
disp(['Interpolating V:',num2str(k),'...']) |
978 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
979 |
Fieldout=inpaint_nans(Field0,0); |
980 |
disp('Done.') |
981 |
|
982 |
disp(['Zeroing V:',num2str(k),' within topography']) |
983 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
984 |
disp('Done.') |
985 |
end |
986 |
|
987 |
for j=1:length(pickout) |
988 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
989 |
Xhere=fout{'X'}(:); |
990 |
Yhere=fout{'Yp1'}(:); |
991 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
992 |
fout{'V'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
993 |
close(fout) |
994 |
end |
995 |
|
996 |
%%%%%%%%%%%%%% gVnm1 |
997 |
|
998 |
Fieldoutk=zeros([length(Zcomp) size(xcout)]); |
999 |
for k=1:length(Zcomp) |
1000 |
Fieldin=NaN*ones(size(xcin)); |
1001 |
Fieldout=NaN*ones(size(xcout)); |
1002 |
|
1003 |
for j=1:length(pickin) |
1004 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
1005 |
gin=netcdf([dirin '/' gridin(j).name],'nowrite'); |
1006 |
Xhere=fin{'X'}(:); |
1007 |
Yhere=fin{'Yp1'}(:); |
1008 |
Fieldinhere=squeeze(fin{'gVnm1'}(snap,k,:,:)); |
1009 |
HFacinhere=squeeze(gin{'HFacS'}(k,:,:)); |
1010 |
for ii=1:length(Xhere) |
1011 |
for jj=1:length(Yhere) |
1012 |
[jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)); |
1013 |
Fieldin(jjj,iii)=Fieldinhere(jj,ii); |
1014 |
if (HFacinhere(jj,ii)==0) |
1015 |
Fieldin(jjj,iii)=NaN; |
1016 |
end |
1017 |
end |
1018 |
end |
1019 |
close(fin) |
1020 |
close(gin) |
1021 |
end |
1022 |
% This utilizes periodic geometry |
1023 |
%Note that there is already a periodic repeat in gVnm1... |
1024 |
Fieldin(:,1)=Fieldin(:,end-1); |
1025 |
Fieldin(:,end)=Fieldin(:,2); |
1026 |
|
1027 |
Fieldin(1,:)=Fieldin(end-2,:); |
1028 |
Fieldin(end,:)=Fieldin(3,:); |
1029 |
|
1030 |
disp(['Interpolating gVnm1:',num2str(k),'...']) |
1031 |
Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
1032 |
Fieldout=inpaint_nans(Field0,0); |
1033 |
disp('Done.') |
1034 |
|
1035 |
disp(['Zeroing gVnm1:',num2str(k),' within topography']) |
1036 |
Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:)); |
1037 |
disp('Done.') |
1038 |
end |
1039 |
|
1040 |
for j=1:length(pickout) |
1041 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
1042 |
Xhere=fout{'X'}(:); |
1043 |
Yhere=fout{'Yp1'}(:); |
1044 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
1045 |
fout{'gVnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
1046 |
close(fout) |
1047 |
end |
1048 |
|
1049 |
|