/[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.4 - (show annotations) (download)
Wed Oct 19 18:47:16 2005 UTC (18 years, 9 months ago) by baylor
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint58e_post, checkpoint58u_post, checkpoint58r_post, checkpoint57y_post, checkpoint58n_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint58j_post, checkpoint58, checkpoint58f_post, checkpoint57x_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint58a_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint58b_post, checkpoint58m_post
Changes since 1.3: +8 -6 lines
Oops.  Convert script back to function.

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

  ViewVC Help
Powered by ViewVC 1.1.22