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

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

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

revision 1.2 by baylor, Mon May 23 20:11:22 2005 UTC revision 1.3 by baylor, Wed Oct 19 18:30:02 2005 UTC
# Line 1  Line 1 
1  function interpickups(dirin,dirout)  %function interpicks(dirin,dirout,varargin)
2  %function interpickups(dirin,dirout)  % function interpicks(dirin,dirout,snap)
3  %This function interpolates the data in  %
4  %a set of mnc pickup files from the MITgcm  % This function interpolates the data in
5  %given in dirin/pickup.*.nc  % a set of mnc pickup files and grid files from the MITgcm
6  %(this can be a global file or a  % given in dirin/pickup.*.nc and in dirin/grid.*.nc
7  % collection of per-processor pickup files)  % (this can be a global file or a collection of per-tile pickup files)
8  %to the pickup files dirout/pickup.*.nc  % to the pickup files dirout/pickup.*.nc and the grid files dirout/grid.*.nc
9  %(this, too, can be a global file or a  % (this, too, can be a global file or a collection of per-tile pickup files)
10  % collection of per-processor pickup files)  %
11  %Extrapolation takes place near the domain edges if domain size  % Extrapolation takes place near the domain edges if domain size
12  %of pickout is larger than that of pickin.  % 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    % 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'])  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'])  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');  fin=netcdf([dirin '/' pickin(1).name],'nowrite');
44    gin=netcdf([dirin '/' gridin(1).name],'nowrite');
45  Zcomp=fin{'Z'}(:);  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'}(:);  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'}(:);  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'}(:);  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'}(:);  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)  close(fin)
77    close(gin)
78    
79  for i=2:length(pickin)  for i=2:length(pickin)
80    fin=netcdf([dirin '/' pickin(i).name],'nowrite');    fin=netcdf([dirin '/' pickin(i).name],'nowrite');
81    Z=fin{'Z'}(:);    Z=fin{'Z'}(:);
82    if sum(Zcomp~=Z)    if (sum(Zcomp~=Z)>0)
83      error('Z','Incompatible vertical axes in input gridfiles')      error('Z','Incompatible vertical axes in input pickups:',num2str(i))
84    end    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'}(:)]);    Xin=sort([Xin;fin{'X'}(:)]);
93    Xp1in=sort([Xp1in;fin{'Xp1'}(:)]);    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'}(:)]);    Yin=sort([Yin;fin{'Y'}(:)]);
103    Yp1in=sort([Yp1in;fin{'Yp1'}(:)]);    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);    close(fin);
113      close(gin);
114  end  end
115    
116  store=[Xin(1)];  store=[Xin(1)];
# Line 43  for i=2:length(Xin) Line 120  for i=2:length(Xin)
120    end    end
121  end  end
122  Xin=store';  Xin=store';
123    clear gXin
124    
125  store=[Xp1in(1)];  store=[Xp1in(1)];
126  for i=2:length(Xp1in)  for i=2:length(Xp1in)
# Line 51  for i=2:length(Xp1in) Line 129  for i=2:length(Xp1in)
129    end    end
130  end  end
131  Xp1in=store';  Xp1in=store';
132    clear gXp1in
133    
134  store=[Yin(1)];  store=[Yin(1)];
135  for i=2:length(Yin)  for i=2:length(Yin)
# Line 59  for i=2:length(Yin) Line 138  for i=2:length(Yin)
138    end    end
139  end  end
140  Yin=store';  Yin=store';
141    clear gYin
142    
143  store=[Yp1in(1)];  store=[Yp1in(1)];
144  for i=2:length(Yp1in)  for i=2:length(Yp1in)
# Line 67  for i=2:length(Yp1in) Line 147  for i=2:length(Yp1in)
147    end    end
148  end  end
149  Yp1in=store';  Yp1in=store';
150    clear gYp1in
151    
152    %%%%%%%%%%%%%%% OUTPUT SANITY
153  fout=netcdf([dirout '/' pickout(1).name],'nowrite');  fout=netcdf([dirout '/' pickout(1).name],'nowrite');
154    gout=netcdf([dirout '/' gridout(1).name],'nowrite');
155  Zcomp=fout{'Z'}(:);  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'}(:);  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'}(:);  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'}(:);  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'}(:);  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)  close(fout)
187    close(gout)
188    
189  for i=2:length(pickout)  for i=2:length(pickout)
190    fout=netcdf([dirout '/' pickout(i).name],'nowrite');    fout=netcdf([dirout '/' pickout(i).name],'nowrite');
191    Z=fout{'Z'}(:);    Z=fout{'Z'}(:);
192    if sum(Zcomp~=Z)    if (sum(Zcomp~=Z)>0)
193      error('Z','Incompatible vertical axes in output gridfiles')      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    end
201    
202    Xout=sort([Xout;fout{'X'}(:)]);    Xout=sort([Xout;fout{'X'}(:)]);
203    Xp1out=sort([Xp1out;fout{'Xp1'}(:)]);    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'}(:)]);    Yout=sort([Yout;fout{'Y'}(:)]);
213    Yp1out=sort([Yp1out;fout{'Yp1'}(:)]);    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);    close(fout);
223      close(gout);
224  end  end
225    
226  store=[Xout(1)];  store=[Xout(1)];
# Line 96  for i=2:length(Xout) Line 230  for i=2:length(Xout)
230    end    end
231  end  end
232  Xout=store';  Xout=store';
233    clear gXout
234    
235  store=[Xp1out(1)];  store=[Xp1out(1)];
236  for i=2:length(Xp1out)  for i=2:length(Xp1out)
# Line 104  for i=2:length(Xp1out) Line 239  for i=2:length(Xp1out)
239    end    end
240  end  end
241  Xp1out=store';  Xp1out=store';
242    clear gXp1out
243    
244  store=[Yout(1)];  store=[Yout(1)];
245  for i=2:length(Yout)  for i=2:length(Yout)
# Line 112  for i=2:length(Yout) Line 248  for i=2:length(Yout)
248    end    end
249  end  end
250  Yout=store';  Yout=store';
251    clear gYout
252    
253  store=[Yp1out(1)];  store=[Yp1out(1)];
254  for i=2:length(Yp1out)  for i=2:length(Yp1out)
# Line 120  for i=2:length(Yp1out) Line 257  for i=2:length(Yp1out)
257    end    end
258  end  end
259  Yp1out=store';  Yp1out=store';
260    clear gYp1out
261    
262    
263  % First, do the Centered variables  % 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)];  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,Yin);  [xcin,ycin]=meshgrid(Xinext,Yinext);
270  [xcout,ycout]=meshgrid(Xout,Yout);  [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));  Fieldin=NaN*ones(size(xcin));
296  Fieldout=NaN*ones(size(xcout));  Fieldout=NaN*ones(size(xcout));
297    
# Line 135  Fieldout=NaN*ones(size(xcout)); Line 299  Fieldout=NaN*ones(size(xcout));
299    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
300    for j=1:length(pickin)    for j=1:length(pickin)
301      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
302        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
303      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
304      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
305      Fieldinhere=fin{'Eta'}(:,:);      Fieldinhere=fin{'Eta'}(:,:);
306        HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
307      for ii=1:length(Xhere)      for ii=1:length(Xhere)
308        for jj=1:length(Yhere)        for jj=1:length(Yhere)
309          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii);          [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        end
315      end      end
316      close(fin)      close(fin)
317        close(gin)
318    end    end
319      % This utilizes channel geometry    % This utilizes periodic geometry
320      Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
321      Fieldin(:,end)=Fieldin(:,2);    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    
     disp(['Interpolating Eta ...'])  
     %Nearests used to eliminate Nans in non-periodic direction  
     Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
     Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});  
     for ii=1:length(Fieldout(:))  
        if isnan(Fieldout(ii))  
          Fieldout(ii)=Fieldoutn(ii);  
        end  
     end  
     disp('Done.')  
335    for j=1:length(pickout)    for j=1:length(pickout)
336      fout=netcdf([dirout '/' pickout(j).name],'write');      fout=netcdf([dirout '/' pickout(j).name],'write');
337      Xhere=fout{'X'}(:);      Xhere=fout{'X'}(:);
# Line 168  Fieldout=NaN*ones(size(xcout)); Line 341  Fieldout=NaN*ones(size(xcout));
341      close(fout)      close(fout)
342    end    end
343    
344    %%%%%%%%%%%%% EtaH
345    
346    Fieldin=NaN*ones(size(xcin));
347    Fieldout=NaN*ones(size(xcout));
348    
349    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
350    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
351    for j=1:length(pickin)    for j=1:length(pickin)
352      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
353        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
354      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
355      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
356      Fieldinhere=fin{'EtaH'}(:,:);      Fieldinhere=fin{'EtaH'}(:,:);
357        HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
358      for ii=1:length(Xhere)      for ii=1:length(Xhere)
359        for jj=1:length(Yhere)        for jj=1:length(Yhere)
360          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii);          [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        end
366      end      end
367      close(fin)      close(fin)
368        close(gin)
369    end    end
370      % This utilizes channel geometry    % This utilizes periodic geometry
371      Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
372      Fieldin(:,end)=Fieldin(:,2);    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    
     disp(['Interpolating EtaH ...'])  
     %Nearests used to eliminate Nans in non-periodic direction  
     Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
     Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});  
     for ii=1:length(Fieldout(:))  
        if isnan(Fieldout(ii))  
          Fieldout(ii)=Fieldoutn(ii);  
        end  
     end  
     disp('Done.')  
386    for j=1:length(pickout)    for j=1:length(pickout)
387      fout=netcdf([dirout '/' pickout(j).name],'write');      fout=netcdf([dirout '/' pickout(j).name],'write');
388      Xhere=fout{'X'}(:);      Xhere=fout{'X'}(:);
# Line 206  Fieldout=NaN*ones(size(xcout)); Line 392  Fieldout=NaN*ones(size(xcout));
392      close(fout)      close(fout)
393    end    end
394    
395    %%%%%%%%%%%%% EtaHdt
396    
397    Fieldin=NaN*ones(size(xcin));
398    Fieldout=NaN*ones(size(xcout));
399    
400    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
401    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
402    for j=1:length(pickin)    for j=1:length(pickin)
403      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
404        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
405      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
406      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
407      Fieldinhere=fin{'dEtaHdt'}(:,:);      Fieldinhere=fin{'dEtaHdt'}(:,:);
408        HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
409      for ii=1:length(Xhere)      for ii=1:length(Xhere)
410        for jj=1:length(Yhere)        for jj=1:length(Yhere)
411          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(jj,ii);          [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        end
417      end      end
418      close(fin)      close(fin)
419        close(gin)
420    end    end
421      % This utilizes channel geometry    % This utilizes periodic geometry
422      Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
423      Fieldin(:,end)=Fieldin(:,2);    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    
     disp(['Interpolating dEtaHdt ...'])  
     %Nearests used to eliminate Nans in non-periodic direction  
     Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
     Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});  
     for ii=1:length(Fieldout(:))  
        if isnan(Fieldout(ii))  
          Fieldout(ii)=Fieldoutn(ii);  
        end  
     end  
     disp('Done.')  
437    for j=1:length(pickout)    for j=1:length(pickout)
438      fout=netcdf([dirout '/' pickout(j).name],'write');      fout=netcdf([dirout '/' pickout(j).name],'write');
439      Xhere=fout{'X'}(:);      Xhere=fout{'X'}(:);
# Line 244  Fieldout=NaN*ones(size(xcout)); Line 443  Fieldout=NaN*ones(size(xcout));
443      close(fout)      close(fout)
444    end    end
445    
446  %S,gSnm1,Temp,gTnm1,phi_nh is on Xc,Rc  %S,gSnm1,Temp,gTnm1,phi_nh are on Xc,Rc
447    
448  Fieldoutk=ones([length(Zcomp) size(xcout)]);  %%%%%%%%%%%%%% S
449  for k=1:length(Z)  
450     Fieldoutk=zeros([length(Zcomp) size(xcout)]);
451     for k=1:length(Zcomp)
452      clear Fieldinhere
453    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
454    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
455    
456    for j=1:length(pickin)    for j=1:length(pickin)
457      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
458        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
459      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
460      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
461      Fieldinhere=fin{'S'}(:,:,:);      Fieldinhere=squeeze(fin{'S'}(snap,k,:,:));
462        HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
463      for ii=1:length(Xhere)      for ii=1:length(Xhere)
464        for jj=1:length(Yhere)        for jj=1:length(Yhere)
465          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
471      end      end
472      close(fin)      close(fin)
473        close(gin)
474    end    end
475    % This utilizes channel geometry    % This utilizes periodic geometry
476    Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
477    Fieldin(:,end)=Fieldin(:,2);    Fieldin(:,end)=Fieldin(:,2);
478    
479    disp(['Interpolating S level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
480    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
481    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
482    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating S:',num2str(k),'...'])
483    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
484       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
485         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
486       end  
487    end    disp(['Zeroing S:',num2str(k),' within topography'])
488    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
489    disp('Done.')    disp('Done.')
490  end  end
491    
492  for j=1:length(pickout)  for j=1:length(pickout)
493    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
494    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);
# Line 287  for j=1:length(pickout) Line 498  for j=1:length(pickout)
498    close(fout)    close(fout)
499  end  end
500    
501    %%%%%%%%%%%%%% gSnm1
502    
503  Fieldoutk=ones([length(Zcomp) size(xcout)]);   Fieldoutk=zeros([length(Zcomp) size(xcout)]);
504  for k=1:length(Z)   for k=1:length(Zcomp)
505    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
506    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
507    
508    for j=1:length(pickin)    for j=1:length(pickin)
509      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
510        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
511      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
512      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
513      Fieldinhere=fin{'gSnm1'}(:,:,:);      Fieldinhere=squeeze(fin{'gSnm1'}(snap,k,:,:));
514        HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
515      for ii=1:length(Xhere)      for ii=1:length(Xhere)
516        for jj=1:length(Yhere)        for jj=1:length(Yhere)
517          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
523      end      end
524      close(fin)      close(fin)
525        close(gin)
526    end    end
527    % This utilizes channel geometry    % This utilizes periodic geometry
528    Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
529    Fieldin(:,end)=Fieldin(:,2);    Fieldin(:,end)=Fieldin(:,2);
530    
531    disp(['Interpolating gSnm1 level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
532    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
533    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
534    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating gSnm1:',num2str(k),'...'])
535    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
536       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
537         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
538       end  
539    end    disp(['Zeroing gSnm1:',num2str(k),' within topography'])
540    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
541    disp('Done.')    disp('Done.')
542  end  end
543    
544  for j=1:length(pickout)  for j=1:length(pickout)
545    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
546    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);
# Line 329  for j=1:length(pickout) Line 550  for j=1:length(pickout)
550    close(fout)    close(fout)
551  end  end
552    
553    %%%%%%%%%%%%%% Temp
554    
555     Fieldoutk=zeros([length(Zcomp) size(xcout)]);
556  Fieldoutk=ones([length(Zcomp) size(xcout)]);   for k=1:length(Zcomp)
 for k=1:length(Z)  
557    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
558    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
559    
560    for j=1:length(pickin)    for j=1:length(pickin)
561      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
562        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
563      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
564      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
565      Fieldinhere=fin{'Temp'}(:,:,:);      Fieldinhere=squeeze(fin{'Temp'}(snap,k,:,:));
566        HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
567      for ii=1:length(Xhere)      for ii=1:length(Xhere)
568        for jj=1:length(Yhere)        for jj=1:length(Yhere)
569          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
575      end      end
576      close(fin)      close(fin)
577        close(gin)
578    end    end
579    % This utilizes channel geometry    % This utilizes periodic geometry
580    Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
581    Fieldin(:,end)=Fieldin(:,2);    Fieldin(:,end)=Fieldin(:,2);
582    
583    disp(['Interpolating Temp level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
584    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
585    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
586    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating Temp:',num2str(k),'...'])
587    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
588       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
589         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
590       end  
591    end    disp(['Zeroing Temp:',num2str(k),' within topography'])
592    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
593    disp('Done.')    disp('Done.')
594  end  end
595    
596  for j=1:length(pickout)  for j=1:length(pickout)
597    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
598    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);
# Line 372  for j=1:length(pickout) Line 602  for j=1:length(pickout)
602    close(fout)    close(fout)
603  end  end
604    
605    %%%%%%%%%%%%%% gTnm1
606    
607  Fieldoutk=ones([length(Zcomp) size(xcout)]);   Fieldoutk=zeros([length(Zcomp) size(xcout)]);
608  for k=1:length(Z)   for k=1:length(Zcomp)
609    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
610    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
611    
612    for j=1:length(pickin)    for j=1:length(pickin)
613      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
614        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
615      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
616      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
617      Fieldinhere=fin{'gTnm1'}(:,:,:);      Fieldinhere=squeeze(fin{'gTnm1'}(snap,k,:,:));
618        HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
619      for ii=1:length(Xhere)      for ii=1:length(Xhere)
620        for jj=1:length(Yhere)        for jj=1:length(Yhere)
621          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
627      end      end
628      close(fin)      close(fin)
629        close(gin)
630    end    end
631    % This utilizes channel geometry    % This utilizes periodic geometry
632    Fieldin(:,1)=Fieldin(:,end-1);    Fieldin(:,1)=Fieldin(:,end-1);
633    Fieldin(:,end)=Fieldin(:,2);    Fieldin(:,end)=Fieldin(:,2);
634    
635    disp(['Interpolating gTnm1 level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
636    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
637    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
638    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating gTnm1:',num2str(k),'...'])
639    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
640       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
641         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
642       end  
643    end    disp(['Zeroing gTnm1:',num2str(k),' within topography'])
644    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
645    disp('Done.')    disp('Done.')
646  end  end
647    
648  for j=1:length(pickout)  for j=1:length(pickout)
649    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
650    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);
# Line 414  for j=1:length(pickout) Line 654  for j=1:length(pickout)
654    close(fout)    close(fout)
655  end  end
656    
657    %%%%%%%%%%%%%% phi_nh
658  fin=netcdf([dirin '/' pickin(1).name],'nowrite');  fin=netcdf([dirin '/' pickin(1).name],'nowrite');
659  status=fin{'phi_nh'};  status=fin{'phi_nh'};
   
660  if ~isempty(status)  if ~isempty(status)
661    Fieldoutk=ones([length(Zcomp) size(xcout)]);  
662    for k=1:length(Z)   Fieldoutk=zeros([length(Zcomp) size(xcout)]);
663      Fieldin=NaN*ones(size(xcin));   for k=1:length(Zcomp)
664      Fieldout=NaN*ones(size(xcout));    Fieldin=NaN*ones(size(xcin));
665      for j=1:length(pickin)    Fieldout=NaN*ones(size(xcout));
666        fin=netcdf([dirin '/' pickin(j).name],'nowrite');  
667        Xhere=fin{'X'}(:);    for j=1:length(pickin)
668        Yhere=fin{'Y'}(:);      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
669        Fieldinhere=fin{'phi_nh'}(:,:,:);      gin=netcdf([dirin '/' gridin(j).name],'nowrite');
670        for ii=1:length(Xhere)      Xhere=fin{'X'}(:);
671          for jj=1:length(Yhere)      Yhere=fin{'Y'}(:);
672            Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);      Fieldinhere=squeeze(fin{'phi_nh'}(snap,k,:,:));
673          end      HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
674        end      for ii=1:length(Xhere)
675        close(fin)        for jj=1:length(Yhere)
676      end          [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
677      % This utilizes channel geometry          Fieldin(jjj,iii)=Fieldinhere(jj,ii);
678      Fieldin(:,1)=Fieldin(:,end-1);          if (HFacinhere(jj,ii)==0)
679      Fieldin(:,end)=Fieldin(:,2);            Fieldin(jjj,iii)=NaN;
680            end
681      disp(['Interpolating phi_nh level=' num2str(k) '...'])        end
     %Nearests used to eliminate Nans in non-periodic direction  
     Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
     Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});  
     for ii=1:length(Fieldout(:))  
        if isnan(Fieldout(ii))  
          Fieldout(ii)=Fieldoutn(ii);  
        end  
682      end      end
683      Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));      close(fin)
684      disp('Done.')      close(gin)
   end  
   for j=1:length(pickout)  
     fout=netcdf([dirout '/' pickout(j).name],'write');  
     Xhere=fout{'X'}(:);  
     Yhere=fout{'Y'}(:);  
     [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));  
     fout{'phi_nh'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));  
     close(fout)  
685    end    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  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');  fin=netcdf([dirin '/' pickin(1).name],'nowrite');
718  status=fin{'gW'};  status=fin{'gW'};
719    
720  if ~isempty(status)  if ~isempty(status)
721    Fieldoutk=ones([length(Zcomp) size(xcout)]);  
722    for k=1:length(Z)   Fieldoutk=zeros([length(Zcomp) size(xcout)]);
723      Fieldin=NaN*ones(size(xcin));   for k=1:length(Zcomp)
724      Fieldout=NaN*ones(size(xcout));    Fieldin=NaN*ones(size(xcin));
725      for j=1:length(pickin)    Fieldout=NaN*ones(size(xcout));
726        fin=netcdf([dirin '/' pickin(j).name],'nowrite');  
727        Xhere=fin{'X'}(:);    for j=1:length(pickin)
728        Yhere=fin{'Y'}(:);      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
729        Fieldinhere=fin{'gW'}(:,:,:);      gin=netcdf([dirin '/' gridin(j).name],'nowrite');
730        for ii=1:length(Xhere)      Xhere=fin{'X'}(:);
731          for jj=1:length(Yhere)      Yhere=fin{'Y'}(:);
732            Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);      Fieldinhere=squeeze(fin{'gW'}(snap,k,:,:));
733          end      HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
734        end      for ii=1:length(Xhere)
735        close(fin)        for jj=1:length(Yhere)
736      end          [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
737      % This utilizes channel geometry          Fieldin(jjj,iii)=Fieldinhere(jj,ii);
738      Fieldin(:,1)=Fieldin(:,end-1);          if (HFacinhere(jj,ii)==0)
739      Fieldin(:,end)=Fieldin(:,2);            Fieldin(jjj,iii)=NaN;
740            end
741      disp(['Interpolating gW level=' num2str(k) '...'])        end
     %Nearests used to eliminate Nans in non-periodic direction  
     Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
     Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});  
     for ii=1:length(Fieldout(:))  
        if isnan(Fieldout(ii))  
          Fieldout(ii)=Fieldoutn(ii);  
        end  
742      end      end
743      Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));      close(fin)
744      disp('Done.')      close(gin)
   end  
   for j=1:length(pickout)  
     fout=netcdf([dirout '/' pickout(j).name],'write');  
     Xhere=fout{'X'}(:);  
     Yhere=fout{'Y'}(:);  
     [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));  
     fout{'gW'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));  
     close(fout)  
745    end    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  end
771    end %isempty(status)
772    
773  %U,gUnm1 is on XU,Rc  %U,gUnm1 is on XU,Rc
774  %Note that there is already a periodic repeat in Uin...  %Note that there is already a periodic repeat in Uin...
775  Xp1inext=[2*Xp1in(1)-Xp1in(2);Xp1in;2*Xp1in(end)-Xp1in(end-1)];  Xp1inext=[2*Xp1in(1)-Xp1in(2);Xp1in;2*Xp1in(end)-Xp1in(end-1)];
776    
777  [xcin,ycin]=meshgrid(Xp1inext,Yin);  [xcin,ycin]=meshgrid(Xp1inext,Yinext);
778  [xcout,ycout]=meshgrid(Xp1out,Yout);  [xcout,ycout]=meshgrid(Xp1out,Yout);
779    
780  Fieldoutk=ones([length(Zcomp) size(xcout)]);  %%%%%%%%%%%%% HFacWoutk
781  for k=1:length(Z)  
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));    Fieldin=NaN*ones(size(xcin));
807    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
808    
809    for j=1:length(pickin)    for j=1:length(pickin)
810      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
811        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
812      Xhere=fin{'Xp1'}(:);      Xhere=fin{'Xp1'}(:);
813      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
814      Fieldinhere=fin{'U'}(:,:,:);      Fieldinhere=squeeze(fin{'U'}(snap,k,:,:));
815        HFacinhere=squeeze(gin{'HFacW'}(k,:,:));
816      for ii=1:length(Xhere)      for ii=1:length(Xhere)
817        for jj=1:length(Yhere)        for jj=1:length(Yhere)
818          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
824      end      end
825      close(fin)      close(fin)
826        close(gin)
827    end    end
828    % This utilizes channel geometry    % This utilizes periodic geometry
829    %Note that there is already a periodic repeat in Uin...    %Note that there is already a periodic repeat in U...
830    Fieldin(:,1)=Fieldin(:,end-2);    Fieldin(:,1)=Fieldin(:,end-2);
831    Fieldin(:,end)=Fieldin(:,3);    Fieldin(:,end)=Fieldin(:,3);
832    
833    disp(['Interpolating U level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
834    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
835    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
836    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating U:',num2str(k),'...'])
837    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
838       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
839         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
840       end  
841    end    disp(['Zeroing U:',num2str(k),' within topography'])
842    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
843    disp('Done.')    disp('Done.')
844  end  end
845    
846  for j=1:length(pickout)  for j=1:length(pickout)
847    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
848    Xhere=fout{'Xp1'}(:);    Xhere=fout{'Xp1'}(:);
# Line 556  for j=1:length(pickout) Line 852  for j=1:length(pickout)
852    close(fout)    close(fout)
853  end  end
854    
855  Fieldoutk=ones([length(Zcomp) size(xcout)]);  %%%%%%%%%%%%%% gUnm1
856  for k=1:length(Z)  
857     Fieldoutk=zeros([length(Zcomp) size(xcout)]);
858     for k=1:length(Zcomp)
859    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
860    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
861    
862    for j=1:length(pickin)    for j=1:length(pickin)
863      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
864        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
865      Xhere=fin{'Xp1'}(:);      Xhere=fin{'Xp1'}(:);
866      Yhere=fin{'Y'}(:);      Yhere=fin{'Y'}(:);
867      Fieldinhere=fin{'gUnm1'}(:,:,:);      Fieldinhere=squeeze(fin{'gUnm1'}(snap,k,:,:));
868        HFacinhere=squeeze(gin{'HFacW'}(k,:,:));
869      for ii=1:length(Xhere)      for ii=1:length(Xhere)
870        for jj=1:length(Yhere)        for jj=1:length(Yhere)
871          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
877      end      end
878      close(fin)      close(fin)
879        close(gin)
880    end    end
881    % This utilizes channel geometry    % This utilizes periodic geometry
882    %Note that there is already a periodic repeat in gUnm1in...    %Note that there is already a periodic repeat in gUnm1...
883    Fieldin(:,1)=Fieldin(:,end-2);    Fieldin(:,1)=Fieldin(:,end-2);
884    Fieldin(:,end)=Fieldin(:,3);    Fieldin(:,end)=Fieldin(:,3);
885    
886    disp(['Interpolating gUnm1 level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-1,:);
887    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(2,:);
888    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
889    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating gUnm1:',num2str(k),'...'])
890    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
891       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
892         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
893       end  
894    end    disp(['Zeroing gUnm1:',num2str(k),' within topography'])
895    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
896    disp('Done.')    disp('Done.')
897  end  end
898    
899  for j=1:length(pickout)  for j=1:length(pickout)
900    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
901    Xhere=fout{'Xp1'}(:);    Xhere=fout{'Xp1'}(:);
# Line 600  end Line 907  end
907    
908  %V,gVnm1 is on XV,Rc  %V,gVnm1 is on XV,Rc
909    
910  [xcin,ycin]=meshgrid(Xinext,Yp1in);  %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);  [xcout,ycout]=meshgrid(Xout,Yp1out);
915    
916  Fieldoutk=ones([length(Zcomp) size(xcout)]);  %%%%%%%%%%%%% HFacSoutk
917  for k=1:length(Z)  
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));    Fieldin=NaN*ones(size(xcin));
943    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
944    
945    for j=1:length(pickin)    for j=1:length(pickin)
946      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
947        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
948      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
949      Yhere=fin{'Yp1'}(:);      Yhere=fin{'Yp1'}(:);
950      Fieldinhere=fin{'V'}(:,:,:);      Fieldinhere=squeeze(fin{'V'}(snap,k,:,:));
951        HFacinhere=squeeze(gin{'HFacS'}(k,:,:));
952      for ii=1:length(Xhere)      for ii=1:length(Xhere)
953        for jj=1:length(Yhere)        for jj=1:length(Yhere)
954          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
960      end      end
961      close(fin)      close(fin)
962        close(gin)
963    end    end
964    % This utilizes channel geometry    % This utilizes periodic geometry
965    %Note that there is already a periodic repeat in Vin...    %Note that there is already a periodic repeat in V...
966    Fieldin(:,1)=Fieldin(:,end-2);    Fieldin(:,1)=Fieldin(:,end-1);
967    Fieldin(:,end)=Fieldin(:,3);    Fieldin(:,end)=Fieldin(:,2);
968    
969    disp(['Interpolating V level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-2,:);
970    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(3,:);
971    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
972    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating V:',num2str(k),'...'])
973    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
974       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
975         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
976       end  
977    end    disp(['Zeroing V:',num2str(k),' within topography'])
978    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
979    disp('Done.')    disp('Done.')
980  end  end
981    
982  for j=1:length(pickout)  for j=1:length(pickout)
983    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
984    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);
# Line 645  for j=1:length(pickout) Line 988  for j=1:length(pickout)
988    close(fout)    close(fout)
989  end  end
990    
991  Fieldoutk=ones([length(Zcomp) size(xcout)]);  %%%%%%%%%%%%%% gVnm1
992  for k=1:length(Z)  
993     Fieldoutk=zeros([length(Zcomp) size(xcout)]);
994     for k=1:length(Zcomp)
995    Fieldin=NaN*ones(size(xcin));    Fieldin=NaN*ones(size(xcin));
996    Fieldout=NaN*ones(size(xcout));    Fieldout=NaN*ones(size(xcout));
997    
998    for j=1:length(pickin)    for j=1:length(pickin)
999      fin=netcdf([dirin '/' pickin(j).name],'nowrite');      fin=netcdf([dirin '/' pickin(j).name],'nowrite');
1000        gin=netcdf([dirin '/' gridin(j).name],'nowrite');
1001      Xhere=fin{'X'}(:);      Xhere=fin{'X'}(:);
1002      Yhere=fin{'Yp1'}(:);      Yhere=fin{'Yp1'}(:);
1003      Fieldinhere=fin{'gVnm1'}(:,:,:);      Fieldinhere=squeeze(fin{'gVnm1'}(snap,k,:,:));
1004        HFacinhere=squeeze(gin{'HFacS'}(k,:,:));
1005      for ii=1:length(Xhere)      for ii=1:length(Xhere)
1006        for jj=1:length(Yhere)        for jj=1:length(Yhere)
1007          Fieldin(find((Xhere(ii)==xcin).*(Yhere(jj)==ycin)))=Fieldinhere(k,jj,ii);          [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        end
1013      end      end
1014      close(fin)      close(fin)
1015        close(gin)
1016    end    end
1017    % This utilizes channel geometry    % This utilizes periodic geometry
1018    %Note that there is already a periodic repeat in gVnm1in...    %Note that there is already a periodic repeat in gVnm1...
1019    Fieldin(:,1)=Fieldin(:,end-2);    Fieldin(:,1)=Fieldin(:,end-1);
1020    Fieldin(:,end)=Fieldin(:,3);    Fieldin(:,end)=Fieldin(:,2);
1021    
1022    disp(['Interpolating gVnm1 level=' num2str(k) '...'])    Fieldin(1,:)=Fieldin(end-2,:);
1023    %Nearests used to eliminate Nans in non-periodic direction    Fieldin(end,:)=Fieldin(3,:);
1024    Fieldoutn=griddata(xcin,ycin,Fieldin,xcout,ycout,'nearest',{'Qt','Qbb','Qc','Qz'});  
1025    Fieldout=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});    disp(['Interpolating gVnm1:',num2str(k),'...'])
1026    for ii=1:length(Fieldout(:))    Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
1027       if isnan(Fieldout(ii))    Fieldout=inpaint_nans(Field0,0);
1028         Fieldout(ii)=Fieldoutn(ii);    disp('Done.')
1029       end  
1030    end    disp(['Zeroing gVnm1:',num2str(k),' within topography'])
1031    Fieldoutk(k,:,:)=reshape(Fieldout,1,length(Fieldout(:,1)),length(Fieldout(1,:)));    Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
1032    disp('Done.')    disp('Done.')
1033  end  end
1034    
1035  for j=1:length(pickout)  for j=1:length(pickout)
1036    fout=netcdf([dirout '/' pickout(j).name],'write');    fout=netcdf([dirout '/' pickout(j).name],'write');
1037    Xhere=fout{'X'}(:);    Xhere=fout{'X'}(:);

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

  ViewVC Help
Powered by ViewVC 1.1.22