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