/[MITgcm]/MITgcm/utils/matlab/interpickups.m
ViewVC logotype

Contents of /MITgcm/utils/matlab/interpickups.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.3 - (show annotations) (download)
Wed Oct 19 18:30:02 2005 UTC (18 years, 9 months ago) by baylor
Branch: MAIN
Changes since 1.2: +624 -270 lines
Adapt interpickups to be much more reliable regarding topography.  Now uses the inpaint_nans.m interpolation (which is now added to repository with the permission of the author) to fill through topography.  Then, new topography is used to zero the fields.  This should allow seamounts or walls to appear and dissapear between pickup files without unexpected consequences.

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

  ViewVC Help
Powered by ViewVC 1.1.22