1 |
function interpickups(dirin,dirout) |
2 |
%function interpickups(dirin,dirout) |
3 |
%This function interpolates the data in |
4 |
%a set of mnc pickup files from the MITgcm |
5 |
%given in dirin/pickup.*.nc |
6 |
%(this can be a global file or a |
7 |
% collection of per-processor pickup files) |
8 |
%to the pickup files dirout/pickup.*.nc |
9 |
%(this, too, can be a global file or a |
10 |
% collection of per-processor pickup files) |
11 |
%Extrapolation takes place near the domain edges if domain size |
12 |
%of pickout is larger than that of pickin. |
13 |
%The number of vertical levels must be the same in the two sets of files. |
14 |
|
15 |
pickin=dir([dirin '/pickup.*.nc']) |
16 |
pickout=dir([dirout '/pickup.*.nc']) |
17 |
|
18 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
19 |
Zcomp=fin{'Z'}(:); |
20 |
Xin=fin{'X'}(:); |
21 |
Yin=fin{'Y'}(:); |
22 |
Xp1in=fin{'Xp1'}(:); |
23 |
Yp1in=fin{'Yp1'}(:); |
24 |
close(fin) |
25 |
|
26 |
for i=2:length(pickin) |
27 |
fin=netcdf([dirin '/' pickin(i).name],'nowrite'); |
28 |
Z=fin{'Z'}(:); |
29 |
if sum(Zcomp~=Z) |
30 |
error('Z','Incompatible vertical axes in input gridfiles') |
31 |
end |
32 |
Xin=sort([Xin;fin{'X'}(:)]); |
33 |
Xp1in=sort([Xp1in;fin{'Xp1'}(:)]); |
34 |
Yin=sort([Yin;fin{'Y'}(:)]); |
35 |
Yp1in=sort([Yp1in;fin{'Yp1'}(:)]); |
36 |
close(fin); |
37 |
end |
38 |
|
39 |
store=[Xin(1)]; |
40 |
for i=2:length(Xin) |
41 |
if Xin(i-1)~=Xin(i) |
42 |
store(end+1)=Xin(i); |
43 |
end |
44 |
end |
45 |
Xin=store'; |
46 |
|
47 |
store=[Xp1in(1)]; |
48 |
for i=2:length(Xp1in) |
49 |
if Xp1in(i-1)~=Xp1in(i) |
50 |
store(end+1)=Xp1in(i); |
51 |
end |
52 |
end |
53 |
Xp1in=store'; |
54 |
|
55 |
store=[Yin(1)]; |
56 |
for i=2:length(Yin) |
57 |
if Yin(i-1)~=Yin(i) |
58 |
store(end+1)=Yin(i); |
59 |
end |
60 |
end |
61 |
Yin=store'; |
62 |
|
63 |
store=[Yp1in(1)]; |
64 |
for i=2:length(Yp1in) |
65 |
if Yp1in(i-1)~=Yp1in(i) |
66 |
store(end+1)=Yp1in(i); |
67 |
end |
68 |
end |
69 |
Yp1in=store'; |
70 |
|
71 |
fout=netcdf([dirout '/' pickout(1).name],'nowrite'); |
72 |
Zcomp=fout{'Z'}(:); |
73 |
Xout=fout{'X'}(:); |
74 |
Yout=fout{'Y'}(:); |
75 |
Xp1out=fout{'Xp1'}(:); |
76 |
Yp1out=fout{'Yp1'}(:); |
77 |
close(fout) |
78 |
|
79 |
for i=2:length(pickout) |
80 |
fout=netcdf([dirout '/' pickout(i).name],'nowrite'); |
81 |
Z=fout{'Z'}(:); |
82 |
if sum(Zcomp~=Z) |
83 |
error('Z','Incompatible vertical axes in output gridfiles') |
84 |
end |
85 |
Xout=sort([Xout;fout{'X'}(:)]); |
86 |
Xp1out=sort([Xp1out;fout{'Xp1'}(:)]); |
87 |
Yout=sort([Yout;fout{'Y'}(:)]); |
88 |
Yp1out=sort([Yp1out;fout{'Yp1'}(:)]); |
89 |
close(fout); |
90 |
end |
91 |
|
92 |
store=[Xout(1)]; |
93 |
for i=2:length(Xout) |
94 |
if Xout(i-1)~=Xout(i) |
95 |
store(end+1)=Xout(i); |
96 |
end |
97 |
end |
98 |
Xout=store'; |
99 |
|
100 |
store=[Xp1out(1)]; |
101 |
for i=2:length(Xp1out) |
102 |
if Xp1out(i-1)~=Xp1out(i) |
103 |
store(end+1)=Xp1out(i); |
104 |
end |
105 |
end |
106 |
Xp1out=store'; |
107 |
|
108 |
store=[Yout(1)]; |
109 |
for i=2:length(Yout) |
110 |
if Yout(i-1)~=Yout(i) |
111 |
store(end+1)=Yout(i); |
112 |
end |
113 |
end |
114 |
Yout=store'; |
115 |
|
116 |
store=[Yp1out(1)]; |
117 |
for i=2:length(Yp1out) |
118 |
if Yp1out(i-1)~=Yp1out(i) |
119 |
store(end+1)=Yp1out(i); |
120 |
end |
121 |
end |
122 |
Yp1out=store'; |
123 |
|
124 |
|
125 |
% First, do the Centered variables |
126 |
Xinext=[2*Xin(1)-Xin(2);Xin;2*Xin(end)-Xin(end-1)]; |
127 |
|
128 |
[xcin,ycin]=meshgrid(Xinext,Yin); |
129 |
[xcout,ycout]=meshgrid(Xout,Yout); |
130 |
|
131 |
Fieldin=NaN*ones(size(xcin)); |
132 |
Fieldout=NaN*ones(size(xcout)); |
133 |
|
134 |
Fieldin=NaN*ones(size(xcin)); |
135 |
Fieldout=NaN*ones(size(xcout)); |
136 |
for j=1:length(pickin) |
137 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
138 |
Xhere=fin{'X'}(:); |
139 |
Yhere=fin{'Y'}(:); |
140 |
Fieldinhere=fin{'Eta'}(:,:); |
141 |
for ii=1:length(Xhere) |
142 |
for jj=1:length(Yhere) |
143 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii); |
144 |
end |
145 |
end |
146 |
close(fin) |
147 |
end |
148 |
% This utilizes channel geometry |
149 |
Fieldin(:,1)=Fieldin(:,end-1); |
150 |
Fieldin(:,end)=Fieldin(:,2); |
151 |
|
152 |
disp(['Interpolating Eta ...']) |
153 |
%Nearests used to eliminate Nans in non-periodic direction |
154 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
155 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
156 |
for ii=1:length(Fieldout(:)) |
157 |
if isnan(Fieldout(ii)) |
158 |
Fieldout(ii)=Fieldoutn(ii); |
159 |
end |
160 |
end |
161 |
disp('Done.') |
162 |
for j=1:length(pickout) |
163 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
164 |
Xhere=fout{'X'}(:); |
165 |
Yhere=fout{'Y'}(:); |
166 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
167 |
fout{'Eta'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
168 |
close(fout) |
169 |
end |
170 |
|
171 |
|
172 |
Fieldin=NaN*ones(size(xcin)); |
173 |
Fieldout=NaN*ones(size(xcout)); |
174 |
for j=1:length(pickin) |
175 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
176 |
Xhere=fin{'X'}(:); |
177 |
Yhere=fin{'Y'}(:); |
178 |
Fieldinhere=fin{'EtaH'}(:,:); |
179 |
for ii=1:length(Xhere) |
180 |
for jj=1:length(Yhere) |
181 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii); |
182 |
end |
183 |
end |
184 |
close(fin) |
185 |
end |
186 |
% This utilizes channel geometry |
187 |
Fieldin(:,1)=Fieldin(:,end-1); |
188 |
Fieldin(:,end)=Fieldin(:,2); |
189 |
|
190 |
disp(['Interpolating EtaH ...']) |
191 |
%Nearests used to eliminate Nans in non-periodic direction |
192 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
193 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
194 |
for ii=1:length(Fieldout(:)) |
195 |
if isnan(Fieldout(ii)) |
196 |
Fieldout(ii)=Fieldoutn(ii); |
197 |
end |
198 |
end |
199 |
disp('Done.') |
200 |
for j=1:length(pickout) |
201 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
202 |
Xhere=fout{'X'}(:); |
203 |
Yhere=fout{'Y'}(:); |
204 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
205 |
fout{'EtaH'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
206 |
close(fout) |
207 |
end |
208 |
|
209 |
|
210 |
Fieldin=NaN*ones(size(xcin)); |
211 |
Fieldout=NaN*ones(size(xcout)); |
212 |
for j=1:length(pickin) |
213 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
214 |
Xhere=fin{'X'}(:); |
215 |
Yhere=fin{'Y'}(:); |
216 |
Fieldinhere=fin{'dEtaHdt'}(:,:); |
217 |
for ii=1:length(Xhere) |
218 |
for jj=1:length(Yhere) |
219 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii); |
220 |
end |
221 |
end |
222 |
close(fin) |
223 |
end |
224 |
% This utilizes channel geometry |
225 |
Fieldin(:,1)=Fieldin(:,end-1); |
226 |
Fieldin(:,end)=Fieldin(:,2); |
227 |
|
228 |
disp(['Interpolating dEtaHdt ...']) |
229 |
%Nearests used to eliminate Nans in non-periodic direction |
230 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
231 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
232 |
for ii=1:length(Fieldout(:)) |
233 |
if isnan(Fieldout(ii)) |
234 |
Fieldout(ii)=Fieldoutn(ii); |
235 |
end |
236 |
end |
237 |
disp('Done.') |
238 |
for j=1:length(pickout) |
239 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
240 |
Xhere=fout{'X'}(:); |
241 |
Yhere=fout{'Y'}(:); |
242 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
243 |
fout{'dEtaHdt'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii)); |
244 |
close(fout) |
245 |
end |
246 |
|
247 |
%S,gSnm1,Temp,gTnm1,phi_nh is on Xc,Rc |
248 |
|
249 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
250 |
for k=1:length(Z) |
251 |
Fieldin=NaN*ones(size(xcin)); |
252 |
Fieldout=NaN*ones(size(xcout)); |
253 |
for j=1:length(pickin) |
254 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
255 |
Xhere=fin{'X'}(:); |
256 |
Yhere=fin{'Y'}(:); |
257 |
Fieldinhere=fin{'S'}(:,:,:); |
258 |
for ii=1:length(Xhere) |
259 |
for jj=1:length(Yhere) |
260 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
261 |
end |
262 |
end |
263 |
close(fin) |
264 |
end |
265 |
% This utilizes channel geometry |
266 |
Fieldin(:,1)=Fieldin(:,end-1); |
267 |
Fieldin(:,end)=Fieldin(:,2); |
268 |
|
269 |
disp(['Interpolating S level=' num2str(k) '...']) |
270 |
%Nearests used to eliminate Nans in non-periodic direction |
271 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
272 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
273 |
for ii=1:length(Fieldout(:)) |
274 |
if isnan(Fieldout(ii)) |
275 |
Fieldout(ii)=Fieldoutn(ii); |
276 |
end |
277 |
end |
278 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
279 |
disp('Done.') |
280 |
end |
281 |
for j=1:length(pickout) |
282 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
283 |
Xhere=fout{'X'}(:); |
284 |
Yhere=fout{'Y'}(:); |
285 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
286 |
fout{'S'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
287 |
close(fout) |
288 |
end |
289 |
|
290 |
|
291 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
292 |
for k=1:length(Z) |
293 |
Fieldin=NaN*ones(size(xcin)); |
294 |
Fieldout=NaN*ones(size(xcout)); |
295 |
for j=1:length(pickin) |
296 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
297 |
Xhere=fin{'X'}(:); |
298 |
Yhere=fin{'Y'}(:); |
299 |
Fieldinhere=fin{'gSnm1'}(:,:,:); |
300 |
for ii=1:length(Xhere) |
301 |
for jj=1:length(Yhere) |
302 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
303 |
end |
304 |
end |
305 |
close(fin) |
306 |
end |
307 |
% This utilizes channel geometry |
308 |
Fieldin(:,1)=Fieldin(:,end-1); |
309 |
Fieldin(:,end)=Fieldin(:,2); |
310 |
|
311 |
disp(['Interpolating gSnm1 level=' num2str(k) '...']) |
312 |
%Nearests used to eliminate Nans in non-periodic direction |
313 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
314 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
315 |
for ii=1:length(Fieldout(:)) |
316 |
if isnan(Fieldout(ii)) |
317 |
Fieldout(ii)=Fieldoutn(ii); |
318 |
end |
319 |
end |
320 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
321 |
disp('Done.') |
322 |
end |
323 |
for j=1:length(pickout) |
324 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
325 |
Xhere=fout{'X'}(:); |
326 |
Yhere=fout{'Y'}(:); |
327 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
328 |
fout{'gSnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
329 |
close(fout) |
330 |
end |
331 |
|
332 |
|
333 |
|
334 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
335 |
for k=1:length(Z) |
336 |
Fieldin=NaN*ones(size(xcin)); |
337 |
Fieldout=NaN*ones(size(xcout)); |
338 |
for j=1:length(pickin) |
339 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
340 |
Xhere=fin{'X'}(:); |
341 |
Yhere=fin{'Y'}(:); |
342 |
Fieldinhere=fin{'Temp'}(:,:,:); |
343 |
for ii=1:length(Xhere) |
344 |
for jj=1:length(Yhere) |
345 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
346 |
end |
347 |
end |
348 |
close(fin) |
349 |
end |
350 |
% This utilizes channel geometry |
351 |
Fieldin(:,1)=Fieldin(:,end-1); |
352 |
Fieldin(:,end)=Fieldin(:,2); |
353 |
|
354 |
disp(['Interpolating Temp level=' num2str(k) '...']) |
355 |
%Nearests used to eliminate Nans in non-periodic direction |
356 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
357 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
358 |
for ii=1:length(Fieldout(:)) |
359 |
if isnan(Fieldout(ii)) |
360 |
Fieldout(ii)=Fieldoutn(ii); |
361 |
end |
362 |
end |
363 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
364 |
disp('Done.') |
365 |
end |
366 |
for j=1:length(pickout) |
367 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
368 |
Xhere=fout{'X'}(:); |
369 |
Yhere=fout{'Y'}(:); |
370 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
371 |
fout{'Temp'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
372 |
close(fout) |
373 |
end |
374 |
|
375 |
|
376 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
377 |
for k=1:length(Z) |
378 |
Fieldin=NaN*ones(size(xcin)); |
379 |
Fieldout=NaN*ones(size(xcout)); |
380 |
for j=1:length(pickin) |
381 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
382 |
Xhere=fin{'X'}(:); |
383 |
Yhere=fin{'Y'}(:); |
384 |
Fieldinhere=fin{'gTnm1'}(:,:,:); |
385 |
for ii=1:length(Xhere) |
386 |
for jj=1:length(Yhere) |
387 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
388 |
end |
389 |
end |
390 |
close(fin) |
391 |
end |
392 |
% This utilizes channel geometry |
393 |
Fieldin(:,1)=Fieldin(:,end-1); |
394 |
Fieldin(:,end)=Fieldin(:,2); |
395 |
|
396 |
disp(['Interpolating gTnm1 level=' num2str(k) '...']) |
397 |
%Nearests used to eliminate Nans in non-periodic direction |
398 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
399 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
400 |
for ii=1:length(Fieldout(:)) |
401 |
if isnan(Fieldout(ii)) |
402 |
Fieldout(ii)=Fieldoutn(ii); |
403 |
end |
404 |
end |
405 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
406 |
disp('Done.') |
407 |
end |
408 |
for j=1:length(pickout) |
409 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
410 |
Xhere=fout{'X'}(:); |
411 |
Yhere=fout{'Y'}(:); |
412 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
413 |
fout{'gTnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
414 |
close(fout) |
415 |
end |
416 |
|
417 |
|
418 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
419 |
status=fin{'phi_nh'}; |
420 |
|
421 |
if ~isempty(status) |
422 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
423 |
for k=1:length(Z) |
424 |
Fieldin=NaN*ones(size(xcin)); |
425 |
Fieldout=NaN*ones(size(xcout)); |
426 |
for j=1:length(pickin) |
427 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
428 |
Xhere=fin{'X'}(:); |
429 |
Yhere=fin{'Y'}(:); |
430 |
Fieldinhere=fin{'phi_nh'}(:,:,:); |
431 |
for ii=1:length(Xhere) |
432 |
for jj=1:length(Yhere) |
433 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
434 |
end |
435 |
end |
436 |
close(fin) |
437 |
end |
438 |
% This utilizes channel geometry |
439 |
Fieldin(:,1)=Fieldin(:,end-1); |
440 |
Fieldin(:,end)=Fieldin(:,2); |
441 |
|
442 |
disp(['Interpolating phi_nh level=' num2str(k) '...']) |
443 |
%Nearests used to eliminate Nans in non-periodic direction |
444 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
445 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
446 |
for ii=1:length(Fieldout(:)) |
447 |
if isnan(Fieldout(ii)) |
448 |
Fieldout(ii)=Fieldoutn(ii); |
449 |
end |
450 |
end |
451 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
452 |
disp('Done.') |
453 |
end |
454 |
for j=1:length(pickout) |
455 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
456 |
Xhere=fout{'X'}(:); |
457 |
Yhere=fout{'Y'}(:); |
458 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
459 |
fout{'phi_nh'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
460 |
close(fout) |
461 |
end |
462 |
end |
463 |
|
464 |
fin=netcdf([dirin '/' pickin(1).name],'nowrite'); |
465 |
status=fin{'gW'}; |
466 |
|
467 |
if ~isempty(status) |
468 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
469 |
for k=1:length(Z) |
470 |
Fieldin=NaN*ones(size(xcin)); |
471 |
Fieldout=NaN*ones(size(xcout)); |
472 |
for j=1:length(pickin) |
473 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
474 |
Xhere=fin{'X'}(:); |
475 |
Yhere=fin{'Y'}(:); |
476 |
Fieldinhere=fin{'gW'}(:,:,:); |
477 |
for ii=1:length(Xhere) |
478 |
for jj=1:length(Yhere) |
479 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
480 |
end |
481 |
end |
482 |
close(fin) |
483 |
end |
484 |
% This utilizes channel geometry |
485 |
Fieldin(:,1)=Fieldin(:,end-1); |
486 |
Fieldin(:,end)=Fieldin(:,2); |
487 |
|
488 |
disp(['Interpolating gW level=' num2str(k) '...']) |
489 |
%Nearests used to eliminate Nans in non-periodic direction |
490 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
491 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
492 |
for ii=1:length(Fieldout(:)) |
493 |
if isnan(Fieldout(ii)) |
494 |
Fieldout(ii)=Fieldoutn(ii); |
495 |
end |
496 |
end |
497 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
498 |
disp('Done.') |
499 |
end |
500 |
for j=1:length(pickout) |
501 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
502 |
Xhere=fout{'X'}(:); |
503 |
Yhere=fout{'Y'}(:); |
504 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
505 |
fout{'gW'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
506 |
close(fout) |
507 |
end |
508 |
end |
509 |
|
510 |
%U,gUnm1 is on XU,Rc |
511 |
%Note that there is already a periodic repeat in Uin... |
512 |
Xp1inext=[2*Xp1in(1)-Xp1in(2);Xp1in;2*Xp1in(end)-Xp1in(end-1)]; |
513 |
|
514 |
[xcin,ycin]=meshgrid(Xp1inext,Yin); |
515 |
[xcout,ycout]=meshgrid(Xp1out,Yout); |
516 |
|
517 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
518 |
for k=1:length(Z) |
519 |
Fieldin=NaN*ones(size(xcin)); |
520 |
Fieldout=NaN*ones(size(xcout)); |
521 |
for j=1:length(pickin) |
522 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
523 |
Xhere=fin{'Xp1'}(:); |
524 |
Yhere=fin{'Y'}(:); |
525 |
Fieldinhere=fin{'U'}(:,:,:); |
526 |
for ii=1:length(Xhere) |
527 |
for jj=1:length(Yhere) |
528 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
529 |
end |
530 |
end |
531 |
close(fin) |
532 |
end |
533 |
% This utilizes channel geometry |
534 |
%Note that there is already a periodic repeat in Uin... |
535 |
Fieldin(:,1)=Fieldin(:,end-2); |
536 |
Fieldin(:,end)=Fieldin(:,3); |
537 |
|
538 |
disp(['Interpolating U level=' num2str(k) '...']) |
539 |
%Nearests used to eliminate Nans in non-periodic direction |
540 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
541 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
542 |
for ii=1:length(Fieldout(:)) |
543 |
if isnan(Fieldout(ii)) |
544 |
Fieldout(ii)=Fieldoutn(ii); |
545 |
end |
546 |
end |
547 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
548 |
disp('Done.') |
549 |
end |
550 |
for j=1:length(pickout) |
551 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
552 |
Xhere=fout{'Xp1'}(:); |
553 |
Yhere=fout{'Y'}(:); |
554 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
555 |
fout{'U'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
556 |
close(fout) |
557 |
end |
558 |
|
559 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
560 |
for k=1:length(Z) |
561 |
Fieldin=NaN*ones(size(xcin)); |
562 |
Fieldout=NaN*ones(size(xcout)); |
563 |
for j=1:length(pickin) |
564 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
565 |
Xhere=fin{'Xp1'}(:); |
566 |
Yhere=fin{'Y'}(:); |
567 |
Fieldinhere=fin{'gUnm1'}(:,:,:); |
568 |
for ii=1:length(Xhere) |
569 |
for jj=1:length(Yhere) |
570 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
571 |
end |
572 |
end |
573 |
close(fin) |
574 |
end |
575 |
% This utilizes channel geometry |
576 |
%Note that there is already a periodic repeat in gUnm1in... |
577 |
Fieldin(:,1)=Fieldin(:,end-2); |
578 |
Fieldin(:,end)=Fieldin(:,3); |
579 |
|
580 |
disp(['Interpolating gUnm1 level=' num2str(k) '...']) |
581 |
%Nearests used to eliminate Nans in non-periodic direction |
582 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
583 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
584 |
for ii=1:length(Fieldout(:)) |
585 |
if isnan(Fieldout(ii)) |
586 |
Fieldout(ii)=Fieldoutn(ii); |
587 |
end |
588 |
end |
589 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
590 |
disp('Done.') |
591 |
end |
592 |
for j=1:length(pickout) |
593 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
594 |
Xhere=fout{'Xp1'}(:); |
595 |
Yhere=fout{'Y'}(:); |
596 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
597 |
fout{'gUnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
598 |
close(fout) |
599 |
end |
600 |
|
601 |
%V,gVnm1 is on XV,Rc |
602 |
|
603 |
[xcin,ycin]=meshgrid(Xinext,Yp1in); |
604 |
[xcout,ycout]=meshgrid(Xout,Yp1out); |
605 |
|
606 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
607 |
for k=1:length(Z) |
608 |
Fieldin=NaN*ones(size(xcin)); |
609 |
Fieldout=NaN*ones(size(xcout)); |
610 |
for j=1:length(pickin) |
611 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
612 |
Xhere=fin{'X'}(:); |
613 |
Yhere=fin{'Yp1'}(:); |
614 |
Fieldinhere=fin{'V'}(:,:,:); |
615 |
for ii=1:length(Xhere) |
616 |
for jj=1:length(Yhere) |
617 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
618 |
end |
619 |
end |
620 |
close(fin) |
621 |
end |
622 |
% This utilizes channel geometry |
623 |
%Note that there is already a periodic repeat in Vin... |
624 |
Fieldin(:,1)=Fieldin(:,end-2); |
625 |
Fieldin(:,end)=Fieldin(:,3); |
626 |
|
627 |
disp(['Interpolating V level=' num2str(k) '...']) |
628 |
%Nearests used to eliminate Nans in non-periodic direction |
629 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
630 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
631 |
for ii=1:length(Fieldout(:)) |
632 |
if isnan(Fieldout(ii)) |
633 |
Fieldout(ii)=Fieldoutn(ii); |
634 |
end |
635 |
end |
636 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
637 |
disp('Done.') |
638 |
end |
639 |
for j=1:length(pickout) |
640 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
641 |
Xhere=fout{'X'}(:); |
642 |
Yhere=fout{'Yp1'}(:); |
643 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
644 |
fout{'V'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
645 |
close(fout) |
646 |
end |
647 |
|
648 |
Fieldoutk=ones([length(Zcomp) size(xcout)]); |
649 |
for k=1:length(Z) |
650 |
Fieldin=NaN*ones(size(xcin)); |
651 |
Fieldout=NaN*ones(size(xcout)); |
652 |
for j=1:length(pickin) |
653 |
fin=netcdf([dirin '/' pickin(j).name],'nowrite'); |
654 |
Xhere=fin{'X'}(:); |
655 |
Yhere=fin{'Yp1'}(:); |
656 |
Fieldinhere=fin{'gVnm1'}(:,:,:); |
657 |
for ii=1:length(Xhere) |
658 |
for jj=1:length(Yhere) |
659 |
Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii); |
660 |
end |
661 |
end |
662 |
close(fin) |
663 |
end |
664 |
% This utilizes channel geometry |
665 |
%Note that there is already a periodic repeat in gVnm1in... |
666 |
Fieldin(:,1)=Fieldin(:,end-2); |
667 |
Fieldin(:,end)=Fieldin(:,3); |
668 |
|
669 |
disp(['Interpolating gVnm1 level=' num2str(k) '...']) |
670 |
%Nearests used to eliminate Nans in non-periodic direction |
671 |
Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'}); |
672 |
Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'}); |
673 |
for ii=1:length(Fieldout(:)) |
674 |
if isnan(Fieldout(ii)) |
675 |
Fieldout(ii)=Fieldoutn(ii); |
676 |
end |
677 |
end |
678 |
Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:))); |
679 |
disp('Done.') |
680 |
end |
681 |
for j=1:length(pickout) |
682 |
fout=netcdf([dirout '/' pickout(j).name],'write'); |
683 |
Xhere=fout{'X'}(:); |
684 |
Yhere=fout{'Yp1'}(:); |
685 |
[jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end))); |
686 |
fout{'gVnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii)); |
687 |
close(fout) |
688 |
end |
689 |
|
690 |
|