Warning, /utils/matlab/interpickups.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
c23374649e Bayl*0001 function interpickups(dirin,dirout,varargin)
0002 % function interpickups(dirin,dirout,snap)
36eeb8dc0c Bayl*0003 %
0004 % This function interpolates the data in
0005 % a set of mnc pickup files and grid files from the MITgcm
0006 % given in dirin/pickup.*.nc and in dirin/grid.*.nc
0007 % (this can be a global file or a collection of per-tile pickup files)
0008 % to the pickup files dirout/pickup.*.nc and the grid files dirout/grid.*.nc
0009 % (this, too, can be a global file or a collection of per-tile pickup files)
0010 %
0011 % Extrapolation takes place near the domain edges if domain size
0012 % of pickout is larger than that of pickin.
0013 %
0014 % The number of vertical levels must be the same in the two sets of files.
0015 %
0016 % Snap is an optional argument if there is more than one timestep
0017 % in the file. The default is 1.
c23374649e Bayl*0018 %
0019 % May be fishy near boundaries if grid is not uniform...
36eeb8dc0c Bayl*0020
c23374649e Bayl*0021 if nargin==2
36eeb8dc0c Bayl*0022 snap=1
c23374649e Bayl*0023 else
0024 snap=varargin{1}
0025 endif
36eeb8dc0c Bayl*0026
0027 if (strcmp(dirin,dirout))
0028 error('dir','You cant use the same input and output directories!')
0029 end
b9c44f1d6b Bayl*0030
0031 pickin=dir([dirin '/pickup.*.nc'])
36eeb8dc0c Bayl*0032 gridin=dir([dirin '/grid.*.nc'])
0033 if length(pickin)~=length(gridin)
0034 error('in','Incompatible number of input pickups and gridfiles')
0035 end
0036
b9c44f1d6b Bayl*0037 pickout=dir([dirout '/pickup.*.nc'])
36eeb8dc0c Bayl*0038 gridout=dir([dirout '/grid.*.nc'])
0039 if length(pickout)~=length(gridout)
0040 error('out','Incompatible number of output pickups and gridfiles')
0041 end
0042
0043 %%%%%%%%%%%%INPUT SANITY
b9c44f1d6b Bayl*0044
0045 fin=netcdf([dirin '/' pickin(1).name],'nowrite');
36eeb8dc0c Bayl*0046 gin=netcdf([dirin '/' gridin(1).name],'nowrite');
b9c44f1d6b Bayl*0047 Zcomp=fin{'Z'}(:);
36eeb8dc0c Bayl*0048 gZcomp=gin{'Z'}(:);
0049 if (sum(Zcomp~=gZcomp)>0)
0050 error('in','Incompatible Z-axis input pickup and gridfile: 1')
0051 end
0052
b9c44f1d6b Bayl*0053 Xin=fin{'X'}(:);
36eeb8dc0c Bayl*0054 gXin=gin{'X'}(:);
0055
0056 if (sum(gXin~=gXin)>0)
0057 error('in','Incompatible x-axis input pickups and gridfile: 1')
0058 end
0059
b9c44f1d6b Bayl*0060 Yin=fin{'Y'}(:);
36eeb8dc0c Bayl*0061 gYin=gin{'Y'}(:);
0062 if (sum(gYin~=gYin)>0)
0063 error('in','Incompatible y-axis input pickups and gridfile: 1')
0064 end
0065
b9c44f1d6b Bayl*0066 Xp1in=fin{'Xp1'}(:);
36eeb8dc0c Bayl*0067 gXp1in=gin{'Xp1'}(:);
0068 if (sum(gXp1in~=gXp1in)>0)
0069 error('in','Incompatible x-axis input pickups and gridfile: 1')
0070 end
0071
b9c44f1d6b Bayl*0072 Yp1in=fin{'Yp1'}(:);
36eeb8dc0c Bayl*0073 gYp1in=gin{'Yp1'}(:);
0074 if (sum(gYp1in~=gYp1in)>0)
0075 error('in','Incompatible y-axis input pickups and gridfile: 1')
0076 end
0077
b9c44f1d6b Bayl*0078 close(fin)
36eeb8dc0c Bayl*0079 close(gin)
b9c44f1d6b Bayl*0080
0081 for i=2:length(pickin)
0082 fin=netcdf([dirin '/' pickin(i).name],'nowrite');
0083 Z=fin{'Z'}(:);
36eeb8dc0c Bayl*0084 if (sum(Zcomp~=Z)>0)
0085 error('Z','Incompatible vertical axes in input pickups:',num2str(i))
0086 end
0087
0088 gin=netcdf([dirin '/' pickin(i).name],'nowrite');
0089 Z=gin{'Z'}(:);
0090 if (sum(Zcomp~=Z)>0)
0091 error('Z','Incompatible vertical axes in input gridfiles:',num2str(i))
b9c44f1d6b Bayl*0092 end
36eeb8dc0c Bayl*0093
b9c44f1d6b Bayl*0094 Xin=sort([Xin;fin{'X'}(:)]);
0095 Xp1in=sort([Xp1in;fin{'Xp1'}(:)]);
36eeb8dc0c Bayl*0096
0097 gXin=sort([gXin;gin{'X'}(:)]);
0098 gXp1in=sort([gXp1in;gin{'Xp1'}(:)]);
0099
0100 if (sum(gXin~=Xin)>0)
0101 error('X','Incompatible x-axes in input files:',num2str(i))
0102 end
0103
b9c44f1d6b Bayl*0104 Yin=sort([Yin;fin{'Y'}(:)]);
0105 Yp1in=sort([Yp1in;fin{'Yp1'}(:)]);
36eeb8dc0c Bayl*0106
0107 gYin=sort([gYin;fin{'Y'}(:)]);
0108 gYp1in=sort([gYp1in;fin{'Yp1'}(:)]);
0109
0110 if (sum(gYin~=Yin)>0)
0111 error('Y','Incompatible y-axes in input files:',num2str(i))
0112 end
0113
b9c44f1d6b Bayl*0114 close(fin);
36eeb8dc0c Bayl*0115 close(gin);
b9c44f1d6b Bayl*0116 end
0117
0118 store=[Xin(1)];
0119 for i=2:length(Xin)
0120 if Xin(i-1)~=Xin(i)
0121 store(end+1)=Xin(i);
0122 end
0123 end
0124 Xin=store';
36eeb8dc0c Bayl*0125 clear gXin
b9c44f1d6b Bayl*0126
0127 store=[Xp1in(1)];
0128 for i=2:length(Xp1in)
0129 if Xp1in(i-1)~=Xp1in(i)
0130 store(end+1)=Xp1in(i);
0131 end
0132 end
0133 Xp1in=store';
36eeb8dc0c Bayl*0134 clear gXp1in
b9c44f1d6b Bayl*0135
0136 store=[Yin(1)];
0137 for i=2:length(Yin)
0138 if Yin(i-1)~=Yin(i)
0139 store(end+1)=Yin(i);
0140 end
0141 end
0142 Yin=store';
36eeb8dc0c Bayl*0143 clear gYin
b9c44f1d6b Bayl*0144
0145 store=[Yp1in(1)];
0146 for i=2:length(Yp1in)
0147 if Yp1in(i-1)~=Yp1in(i)
0148 store(end+1)=Yp1in(i);
0149 end
0150 end
0151 Yp1in=store';
36eeb8dc0c Bayl*0152 clear gYp1in
b9c44f1d6b Bayl*0153
36eeb8dc0c Bayl*0154 %%%%%%%%%%%%%%% OUTPUT SANITY
b9c44f1d6b Bayl*0155 fout=netcdf([dirout '/' pickout(1).name],'nowrite');
36eeb8dc0c Bayl*0156 gout=netcdf([dirout '/' gridout(1).name],'nowrite');
b9c44f1d6b Bayl*0157 Zcomp=fout{'Z'}(:);
36eeb8dc0c Bayl*0158 gZcomp=gout{'Z'}(:);
0159 if (sum(Zcomp~=gZcomp)>0)
0160 error('out','Incompatible Z-axis output pickup and gridfile: 1')
0161 end
0162
b9c44f1d6b Bayl*0163 Xout=fout{'X'}(:);
36eeb8dc0c Bayl*0164 gXout=gout{'X'}(:);
0165
0166 if (sum(gXout~=gXout)>0)
0167 error('out','Incompatible x-axis output pickups and gridfile: 1')
0168 end
0169
b9c44f1d6b Bayl*0170 Yout=fout{'Y'}(:);
36eeb8dc0c Bayl*0171 gYout=gout{'Y'}(:);
0172 if (sum(gYout~=gYout)>0)
0173 error('out','Incompatible y-axis output pickups and gridfile: 1')
0174 end
0175
b9c44f1d6b Bayl*0176 Xp1out=fout{'Xp1'}(:);
36eeb8dc0c Bayl*0177 gXp1out=gout{'Xp1'}(:);
0178 if (sum(gXp1out~=gXp1out)>0)
0179 error('out','Incompatible x-axis output pickups and gridfile: 1')
0180 end
0181
b9c44f1d6b Bayl*0182 Yp1out=fout{'Yp1'}(:);
36eeb8dc0c Bayl*0183 gYp1out=gout{'Yp1'}(:);
0184 if (sum(gYp1out~=gYp1out)>0)
0185 error('out','Incompatible y-axis output pickups and gridfile: 1')
0186 end
0187
b9c44f1d6b Bayl*0188 close(fout)
36eeb8dc0c Bayl*0189 close(gout)
b9c44f1d6b Bayl*0190
0191 for i=2:length(pickout)
0192 fout=netcdf([dirout '/' pickout(i).name],'nowrite');
0193 Z=fout{'Z'}(:);
36eeb8dc0c Bayl*0194 if (sum(Zcomp~=Z)>0)
0195 error('Z','Incompatible vertical axes in output pickups:',num2str(i))
0196 end
0197
0198 gout=netcdf([dirout '/' pickout(i).name],'nowrite');
0199 Z=gout{'Z'}(:);
0200 if (sum(Zcomp~=Z)>0)
0201 error('Z','Incompatible vertical axes in output gridfiles:',num2str(i))
b9c44f1d6b Bayl*0202 end
36eeb8dc0c Bayl*0203
b9c44f1d6b Bayl*0204 Xout=sort([Xout;fout{'X'}(:)]);
0205 Xp1out=sort([Xp1out;fout{'Xp1'}(:)]);
36eeb8dc0c Bayl*0206
0207 gXout=sort([gXout;gout{'X'}(:)]);
0208 gXp1out=sort([gXp1out;gout{'Xp1'}(:)]);
0209
0210 if (sum(gXout~=Xout)>0)
0211 error('X','Incompatible x-axes in output files:',num2str(i))
0212 end
0213
b9c44f1d6b Bayl*0214 Yout=sort([Yout;fout{'Y'}(:)]);
0215 Yp1out=sort([Yp1out;fout{'Yp1'}(:)]);
36eeb8dc0c Bayl*0216
0217 gYout=sort([gYout;fout{'Y'}(:)]);
0218 gYp1out=sort([gYp1out;fout{'Yp1'}(:)]);
0219
0220 if (sum(gYout~=Yout)>0)
0221 error('Y','Incompatible y-axes in output files:',num2str(i))
0222 end
0223
b9c44f1d6b Bayl*0224 close(fout);
36eeb8dc0c Bayl*0225 close(gout);
b9c44f1d6b Bayl*0226 end
0227
0228 store=[Xout(1)];
0229 for i=2:length(Xout)
0230 if Xout(i-1)~=Xout(i)
0231 store(end+1)=Xout(i);
0232 end
0233 end
0234 Xout=store';
36eeb8dc0c Bayl*0235 clear gXout
b9c44f1d6b Bayl*0236
0237 store=[Xp1out(1)];
0238 for i=2:length(Xp1out)
0239 if Xp1out(i-1)~=Xp1out(i)
0240 store(end+1)=Xp1out(i);
0241 end
0242 end
0243 Xp1out=store';
36eeb8dc0c Bayl*0244 clear gXp1out
b9c44f1d6b Bayl*0245
0246 store=[Yout(1)];
0247 for i=2:length(Yout)
0248 if Yout(i-1)~=Yout(i)
0249 store(end+1)=Yout(i);
0250 end
0251 end
0252 Yout=store';
36eeb8dc0c Bayl*0253 clear gYout
b9c44f1d6b Bayl*0254
0255 store=[Yp1out(1)];
0256 for i=2:length(Yp1out)
0257 if Yp1out(i-1)~=Yp1out(i)
0258 store(end+1)=Yp1out(i);
0259 end
0260 end
0261 Yp1out=store';
36eeb8dc0c Bayl*0262 clear gYp1out
b9c44f1d6b Bayl*0263
0264
0265 % First, do the Centered variables
36eeb8dc0c Bayl*0266
0267 % We assume periodicity as MITgcm does
b9c44f1d6b Bayl*0268 Xinext=[2*Xin(1)-Xin(2);Xin;2*Xin(end)-Xin(end-1)];
36eeb8dc0c Bayl*0269 Yinext=[2*Yin(1)-Yin(2);Yin;2*Yin(end)-Yin(end-1)];
b9c44f1d6b Bayl*0270
36eeb8dc0c Bayl*0271 [xcin,ycin]=meshgrid(Xinext,Yinext);
b9c44f1d6b Bayl*0272 [xcout,ycout]=meshgrid(Xout,Yout);
0273
36eeb8dc0c Bayl*0274 %%%%%%%%%%%%% HFacCoutk
0275
0276 HFacoutk=ones([length(Zcomp) size(xcout)]);
0277
0278 disp(['Calculating HFacC...'])
0279 for j=1:length(gridout)
0280 gout=netcdf([dirout '/' gridout(j).name],'nowrite');
0281 Xhere=gout{'X'}(:);
0282 Yhere=gout{'Y'}(:);
0283 HFacouthere=gout{'HFacC'}(:,:,:);
0284 for ii=1:length(Xhere)
0285 for jj=1:length(Yhere)
0286 [jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout));
0287 HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii);
0288 end
0289 end
0290 close(gout)
0291 end
0292 clear HFacouthere
0293 disp(['Done.'])
0294
0295 %%%%%%%%%%%%% ETA
0296
b9c44f1d6b Bayl*0297 Fieldin=NaN*ones(size(xcin));
0298 Fieldout=NaN*ones(size(xcout));
0299
0300 Fieldin=NaN*ones(size(xcin));
0301 Fieldout=NaN*ones(size(xcout));
0302 for j=1:length(pickin)
0303 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0304 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0305 Xhere=fin{'X'}(:);
0306 Yhere=fin{'Y'}(:);
0307 Fieldinhere=fin{'Eta'}(:,:);
36eeb8dc0c Bayl*0308 HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
b9c44f1d6b Bayl*0309 for ii=1:length(Xhere)
0310 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0311 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0312 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0313 if (HFacinhere(jj,ii)==0)
0314 Fieldin(jjj,iii)=NaN;
0315 end
b9c44f1d6b Bayl*0316 end
0317 end
0318 close(fin)
36eeb8dc0c Bayl*0319 close(gin)
b9c44f1d6b Bayl*0320 end
36eeb8dc0c Bayl*0321 % This utilizes periodic geometry
0322 Fieldin(:,1)=Fieldin(:,end-1);
0323 Fieldin(:,end)=Fieldin(:,2);
0324
0325 Fieldin(1,:)=Fieldin(end-1,:);
0326 Fieldin(end,:)=Fieldin(2,:);
0327
0328 disp(['Interpolating Eta ...'])
0329 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0330 Fieldout=inpaint_nans(Field0,0);
0331 disp('Done.')
0332
0333 disp(['Zeroing Eta within topography'])
0334 Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:));
0335 disp('Done.')
0336
b9c44f1d6b Bayl*0337 for j=1:length(pickout)
0338 fout=netcdf([dirout '/' pickout(j).name],'write');
0339 Xhere=fout{'X'}(:);
0340 Yhere=fout{'Y'}(:);
0341 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0342 fout{'Eta'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii));
0343 close(fout)
0344 end
0345
36eeb8dc0c Bayl*0346 %%%%%%%%%%%%% EtaH
0347
0348 Fieldin=NaN*ones(size(xcin));
0349 Fieldout=NaN*ones(size(xcout));
b9c44f1d6b Bayl*0350
0351 Fieldin=NaN*ones(size(xcin));
0352 Fieldout=NaN*ones(size(xcout));
0353 for j=1:length(pickin)
0354 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0355 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0356 Xhere=fin{'X'}(:);
0357 Yhere=fin{'Y'}(:);
0358 Fieldinhere=fin{'EtaH'}(:,:);
36eeb8dc0c Bayl*0359 HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
b9c44f1d6b Bayl*0360 for ii=1:length(Xhere)
0361 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0362 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0363 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0364 if (HFacinhere(jj,ii)==0)
0365 Fieldin(jjj,iii)=NaN;
0366 end
b9c44f1d6b Bayl*0367 end
0368 end
0369 close(fin)
36eeb8dc0c Bayl*0370 close(gin)
b9c44f1d6b Bayl*0371 end
36eeb8dc0c Bayl*0372 % This utilizes periodic geometry
0373 Fieldin(:,1)=Fieldin(:,end-1);
0374 Fieldin(:,end)=Fieldin(:,2);
0375
0376 Fieldin(1,:)=Fieldin(end-1,:);
0377 Fieldin(end,:)=Fieldin(2,:);
0378
0379 disp(['Interpolating EtaH ...'])
0380 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0381 Fieldout=inpaint_nans(Field0,0);
0382 disp('Done.')
0383
0384 disp(['Zeroing EtaH within topography'])
0385 Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:));
0386 disp('Done.')
0387
b9c44f1d6b Bayl*0388 for j=1:length(pickout)
0389 fout=netcdf([dirout '/' pickout(j).name],'write');
0390 Xhere=fout{'X'}(:);
0391 Yhere=fout{'Y'}(:);
0392 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0393 fout{'EtaH'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii));
0394 close(fout)
0395 end
0396
36eeb8dc0c Bayl*0397 %%%%%%%%%%%%% EtaHdt
0398
0399 Fieldin=NaN*ones(size(xcin));
0400 Fieldout=NaN*ones(size(xcout));
b9c44f1d6b Bayl*0401
0402 Fieldin=NaN*ones(size(xcin));
0403 Fieldout=NaN*ones(size(xcout));
0404 for j=1:length(pickin)
0405 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0406 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0407 Xhere=fin{'X'}(:);
0408 Yhere=fin{'Y'}(:);
0409 Fieldinhere=fin{'dEtaHdt'}(:,:);
36eeb8dc0c Bayl*0410 HFacinhere=squeeze(gin{'HFacC'}(1,:,:));
b9c44f1d6b Bayl*0411 for ii=1:length(Xhere)
0412 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0413 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0414 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0415 if (HFacinhere(jj,ii)==0)
0416 Fieldin(jjj,iii)=NaN;
0417 end
b9c44f1d6b Bayl*0418 end
0419 end
0420 close(fin)
36eeb8dc0c Bayl*0421 close(gin)
b9c44f1d6b Bayl*0422 end
36eeb8dc0c Bayl*0423 % This utilizes periodic geometry
0424 Fieldin(:,1)=Fieldin(:,end-1);
0425 Fieldin(:,end)=Fieldin(:,2);
0426
0427 Fieldin(1,:)=Fieldin(end-1,:);
0428 Fieldin(end,:)=Fieldin(2,:);
0429
0430 disp(['Interpolating dEtaHdt ...'])
0431 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0432 Fieldout=inpaint_nans(Field0,0);
0433 disp('Done.')
0434
0435 disp(['Zeroing dEtaHdt within topography'])
0436 Fieldout=Fieldout.*squeeze(HFacoutk(1,:,:));
0437 disp('Done.')
0438
b9c44f1d6b Bayl*0439 for j=1:length(pickout)
0440 fout=netcdf([dirout '/' pickout(j).name],'write');
0441 Xhere=fout{'X'}(:);
0442 Yhere=fout{'Y'}(:);
0443 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0444 fout{'dEtaHdt'}(:,:)=Fieldout(min(jj):max(jj),min(ii):max(ii));
0445 close(fout)
0446 end
0447
36eeb8dc0c Bayl*0448 %S,gSnm1,Temp,gTnm1,phi_nh are on Xc,Rc
b9c44f1d6b Bayl*0449
36eeb8dc0c Bayl*0450 %%%%%%%%%%%%%% S
0451
0452 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0453 for k=1:length(Zcomp)
0454 clear Fieldinhere
b9c44f1d6b Bayl*0455 Fieldin=NaN*ones(size(xcin));
0456 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0457
b9c44f1d6b Bayl*0458 for j=1:length(pickin)
0459 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0460 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0461 Xhere=fin{'X'}(:);
0462 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0463 Fieldinhere=squeeze(fin{'S'}(snap,k,:,:));
0464 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
b9c44f1d6b Bayl*0465 for ii=1:length(Xhere)
0466 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0467 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0468 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0469 if (HFacinhere(jj,ii)==0)
0470 Fieldin(jjj,iii)=NaN;
0471 end
b9c44f1d6b Bayl*0472 end
0473 end
0474 close(fin)
36eeb8dc0c Bayl*0475 close(gin)
b9c44f1d6b Bayl*0476 end
36eeb8dc0c Bayl*0477 % This utilizes periodic geometry
b9c44f1d6b Bayl*0478 Fieldin(:,1)=Fieldin(:,end-1);
0479 Fieldin(:,end)=Fieldin(:,2);
0480
36eeb8dc0c Bayl*0481 Fieldin(1,:)=Fieldin(end-1,:);
0482 Fieldin(end,:)=Fieldin(2,:);
0483
0484 disp(['Interpolating S:',num2str(k),'...'])
0485 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0486 Fieldout=inpaint_nans(Field0,0);
0487 disp('Done.')
0488
0489 disp(['Zeroing S:',num2str(k),' within topography'])
0490 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0491 disp('Done.')
0492 end
36eeb8dc0c Bayl*0493
b9c44f1d6b Bayl*0494 for j=1:length(pickout)
0495 fout=netcdf([dirout '/' pickout(j).name],'write');
0496 Xhere=fout{'X'}(:);
0497 Yhere=fout{'Y'}(:);
0498 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0499 fout{'S'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0500 close(fout)
0501 end
0502
36eeb8dc0c Bayl*0503 %%%%%%%%%%%%%% gSnm1
b9c44f1d6b Bayl*0504
36eeb8dc0c Bayl*0505 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0506 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0507 Fieldin=NaN*ones(size(xcin));
0508 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0509
b9c44f1d6b Bayl*0510 for j=1:length(pickin)
0511 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0512 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0513 Xhere=fin{'X'}(:);
0514 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0515 Fieldinhere=squeeze(fin{'gSnm1'}(snap,k,:,:));
0516 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
b9c44f1d6b Bayl*0517 for ii=1:length(Xhere)
0518 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0519 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0520 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0521 if (HFacinhere(jj,ii)==0)
0522 Fieldin(jjj,iii)=NaN;
0523 end
b9c44f1d6b Bayl*0524 end
0525 end
0526 close(fin)
36eeb8dc0c Bayl*0527 close(gin)
b9c44f1d6b Bayl*0528 end
36eeb8dc0c Bayl*0529 % This utilizes periodic geometry
b9c44f1d6b Bayl*0530 Fieldin(:,1)=Fieldin(:,end-1);
0531 Fieldin(:,end)=Fieldin(:,2);
0532
36eeb8dc0c Bayl*0533 Fieldin(1,:)=Fieldin(end-1,:);
0534 Fieldin(end,:)=Fieldin(2,:);
0535
0536 disp(['Interpolating gSnm1:',num2str(k),'...'])
0537 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0538 Fieldout=inpaint_nans(Field0,0);
0539 disp('Done.')
0540
0541 disp(['Zeroing gSnm1:',num2str(k),' within topography'])
0542 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0543 disp('Done.')
0544 end
36eeb8dc0c Bayl*0545
b9c44f1d6b Bayl*0546 for j=1:length(pickout)
0547 fout=netcdf([dirout '/' pickout(j).name],'write');
0548 Xhere=fout{'X'}(:);
0549 Yhere=fout{'Y'}(:);
0550 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0551 fout{'gSnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0552 close(fout)
0553 end
0554
36eeb8dc0c Bayl*0555 %%%%%%%%%%%%%% Temp
b9c44f1d6b Bayl*0556
36eeb8dc0c Bayl*0557 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0558 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0559 Fieldin=NaN*ones(size(xcin));
0560 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0561
b9c44f1d6b Bayl*0562 for j=1:length(pickin)
0563 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0564 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0565 Xhere=fin{'X'}(:);
0566 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0567 Fieldinhere=squeeze(fin{'Temp'}(snap,k,:,:));
0568 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
b9c44f1d6b Bayl*0569 for ii=1:length(Xhere)
0570 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0571 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0572 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0573 if (HFacinhere(jj,ii)==0)
0574 Fieldin(jjj,iii)=NaN;
0575 end
b9c44f1d6b Bayl*0576 end
0577 end
0578 close(fin)
36eeb8dc0c Bayl*0579 close(gin)
b9c44f1d6b Bayl*0580 end
36eeb8dc0c Bayl*0581 % This utilizes periodic geometry
b9c44f1d6b Bayl*0582 Fieldin(:,1)=Fieldin(:,end-1);
0583 Fieldin(:,end)=Fieldin(:,2);
0584
36eeb8dc0c Bayl*0585 Fieldin(1,:)=Fieldin(end-1,:);
0586 Fieldin(end,:)=Fieldin(2,:);
0587
0588 disp(['Interpolating Temp:',num2str(k),'...'])
0589 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0590 Fieldout=inpaint_nans(Field0,0);
0591 disp('Done.')
0592
0593 disp(['Zeroing Temp:',num2str(k),' within topography'])
0594 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0595 disp('Done.')
0596 end
36eeb8dc0c Bayl*0597
b9c44f1d6b Bayl*0598 for j=1:length(pickout)
0599 fout=netcdf([dirout '/' pickout(j).name],'write');
0600 Xhere=fout{'X'}(:);
0601 Yhere=fout{'Y'}(:);
0602 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0603 fout{'Temp'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0604 close(fout)
0605 end
0606
36eeb8dc0c Bayl*0607 %%%%%%%%%%%%%% gTnm1
b9c44f1d6b Bayl*0608
36eeb8dc0c Bayl*0609 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0610 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0611 Fieldin=NaN*ones(size(xcin));
0612 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0613
b9c44f1d6b Bayl*0614 for j=1:length(pickin)
0615 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0616 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0617 Xhere=fin{'X'}(:);
0618 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0619 Fieldinhere=squeeze(fin{'gTnm1'}(snap,k,:,:));
0620 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
b9c44f1d6b Bayl*0621 for ii=1:length(Xhere)
0622 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0623 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0624 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0625 if (HFacinhere(jj,ii)==0)
0626 Fieldin(jjj,iii)=NaN;
0627 end
b9c44f1d6b Bayl*0628 end
0629 end
0630 close(fin)
36eeb8dc0c Bayl*0631 close(gin)
b9c44f1d6b Bayl*0632 end
36eeb8dc0c Bayl*0633 % This utilizes periodic geometry
b9c44f1d6b Bayl*0634 Fieldin(:,1)=Fieldin(:,end-1);
0635 Fieldin(:,end)=Fieldin(:,2);
0636
36eeb8dc0c Bayl*0637 Fieldin(1,:)=Fieldin(end-1,:);
0638 Fieldin(end,:)=Fieldin(2,:);
0639
0640 disp(['Interpolating gTnm1:',num2str(k),'...'])
0641 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0642 Fieldout=inpaint_nans(Field0,0);
0643 disp('Done.')
0644
0645 disp(['Zeroing gTnm1:',num2str(k),' within topography'])
0646 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0647 disp('Done.')
0648 end
36eeb8dc0c Bayl*0649
b9c44f1d6b Bayl*0650 for j=1:length(pickout)
0651 fout=netcdf([dirout '/' pickout(j).name],'write');
0652 Xhere=fout{'X'}(:);
0653 Yhere=fout{'Y'}(:);
0654 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0655 fout{'gTnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0656 close(fout)
0657 end
0658
36eeb8dc0c Bayl*0659 %%%%%%%%%%%%%% phi_nh
b9c44f1d6b Bayl*0660 fin=netcdf([dirin '/' pickin(1).name],'nowrite');
0661 status=fin{'phi_nh'};
0662 if ~isempty(status)
36eeb8dc0c Bayl*0663
0664 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0665 for k=1:length(Zcomp)
0666 Fieldin=NaN*ones(size(xcin));
0667 Fieldout=NaN*ones(size(xcout));
0668
0669 for j=1:length(pickin)
0670 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
0671 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
0672 Xhere=fin{'X'}(:);
0673 Yhere=fin{'Y'}(:);
0674 Fieldinhere=squeeze(fin{'phi_nh'}(snap,k,:,:));
0675 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
0676 for ii=1:length(Xhere)
0677 for jj=1:length(Yhere)
0678 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0679 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0680 if (HFacinhere(jj,ii)==0)
0681 Fieldin(jjj,iii)=NaN;
b9c44f1d6b Bayl*0682 end
0683 end
0684 end
36eeb8dc0c Bayl*0685 close(fin)
0686 close(gin)
b9c44f1d6b Bayl*0687 end
36eeb8dc0c Bayl*0688 % This utilizes periodic geometry
0689 Fieldin(:,1)=Fieldin(:,end-1);
0690 Fieldin(:,end)=Fieldin(:,2);
0691
0692 Fieldin(1,:)=Fieldin(end-1,:);
0693 Fieldin(end,:)=Fieldin(2,:);
0694
0695 disp(['Interpolating phi_nh:',num2str(k),'...'])
0696 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0697 Fieldout=inpaint_nans(Field0,0);
0698 disp('Done.')
0699
0700 disp(['Zeroing phi_nh:',num2str(k),' within topography'])
0701 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
0702 disp('Done.')
0703 end
0704
0705 for j=1:length(pickout)
0706 fout=netcdf([dirout '/' pickout(j).name],'write');
0707 Xhere=fout{'X'}(:);
0708 Yhere=fout{'Y'}(:);
0709 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0710 fout{'phi_nh'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0711 close(fout)
b9c44f1d6b Bayl*0712 end
0713
36eeb8dc0c Bayl*0714 end %isempty(status)
0715
0716 %%%%%%%%%%%%%%%%%% gW
0717 %BFK Not sure this is right, since HFacC isn't on the right grid for gW
0718
b9c44f1d6b Bayl*0719 fin=netcdf([dirin '/' pickin(1).name],'nowrite');
0720 status=fin{'gW'};
0721
0722 if ~isempty(status)
36eeb8dc0c Bayl*0723
0724 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0725 for k=1:length(Zcomp)
0726 Fieldin=NaN*ones(size(xcin));
0727 Fieldout=NaN*ones(size(xcout));
0728
0729 for j=1:length(pickin)
0730 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
0731 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
0732 Xhere=fin{'X'}(:);
0733 Yhere=fin{'Y'}(:);
0734 Fieldinhere=squeeze(fin{'gW'}(snap,k,:,:));
0735 HFacinhere=squeeze(gin{'HFacC'}(k,:,:));
0736 for ii=1:length(Xhere)
0737 for jj=1:length(Yhere)
0738 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0739 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0740 if (HFacinhere(jj,ii)==0)
0741 Fieldin(jjj,iii)=NaN;
b9c44f1d6b Bayl*0742 end
0743 end
0744 end
36eeb8dc0c Bayl*0745 close(fin)
0746 close(gin)
b9c44f1d6b Bayl*0747 end
36eeb8dc0c Bayl*0748 % This utilizes periodic geometry
0749 Fieldin(:,1)=Fieldin(:,end-1);
0750 Fieldin(:,end)=Fieldin(:,2);
0751
0752 Fieldin(1,:)=Fieldin(end-1,:);
0753 Fieldin(end,:)=Fieldin(2,:);
0754
0755 disp(['Interpolating gW:',num2str(k),'...'])
0756 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0757 Fieldout=inpaint_nans(Field0,0);
0758 disp('Done.')
0759
0760 disp(['Zeroing gW:',num2str(k),' within topography'])
0761 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
0762 disp('Done.')
b9c44f1d6b Bayl*0763 end
0764
36eeb8dc0c Bayl*0765 for j=1:length(pickout)
0766 fout=netcdf([dirout '/' pickout(j).name],'write');
0767 Xhere=fout{'X'}(:);
0768 Yhere=fout{'Y'}(:);
0769 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0770 fout{'gW'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0771 close(fout)
0772 end
0773 end %isempty(status)
0774
b9c44f1d6b Bayl*0775 %U,gUnm1 is on XU,Rc
0776 %Note that there is already a periodic repeat in Uin...
0777 Xp1inext=[2*Xp1in(1)-Xp1in(2);Xp1in;2*Xp1in(end)-Xp1in(end-1)];
0778
36eeb8dc0c Bayl*0779 [xcin,ycin]=meshgrid(Xp1inext,Yinext);
b9c44f1d6b Bayl*0780 [xcout,ycout]=meshgrid(Xp1out,Yout);
0781
36eeb8dc0c Bayl*0782 %%%%%%%%%%%%% HFacWoutk
0783
0784 HFacout=ones(size(xcout));
0785 HFacoutk=ones([length(Zcomp) size(xcout)]);
0786
0787 disp(['Calculating HFacW...'])
0788 for j=1:length(gridout)
0789 gout=netcdf([dirout '/' gridout(j).name],'nowrite');
0790 Xhere=gout{'Xp1'}(:);
0791 Yhere=gout{'Y'}(:);
0792 HFacouthere=gout{'HFacW'}(:,:,:);
0793 for ii=1:length(Xhere)
0794 for jj=1:length(Yhere)
0795 [jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout));
0796 HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii);
0797 end
0798 end
0799 close(gout)
0800 end
0801 clear HFacouthere
0802 disp(['Done.'])
0803
0804 %%%%%%%%%%%%%% U
0805
0806 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0807 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0808 Fieldin=NaN*ones(size(xcin));
0809 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0810
b9c44f1d6b Bayl*0811 for j=1:length(pickin)
0812 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0813 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0814 Xhere=fin{'Xp1'}(:);
0815 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0816 Fieldinhere=squeeze(fin{'U'}(snap,k,:,:));
0817 HFacinhere=squeeze(gin{'HFacW'}(k,:,:));
b9c44f1d6b Bayl*0818 for ii=1:length(Xhere)
0819 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0820 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0821 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0822 if (HFacinhere(jj,ii)==0)
0823 Fieldin(jjj,iii)=NaN;
0824 end
b9c44f1d6b Bayl*0825 end
0826 end
0827 close(fin)
36eeb8dc0c Bayl*0828 close(gin)
b9c44f1d6b Bayl*0829 end
36eeb8dc0c Bayl*0830 % This utilizes periodic geometry
0831 %Note that there is already a periodic repeat in U...
b9c44f1d6b Bayl*0832 Fieldin(:,1)=Fieldin(:,end-2);
0833 Fieldin(:,end)=Fieldin(:,3);
0834
36eeb8dc0c Bayl*0835 Fieldin(1,:)=Fieldin(end-1,:);
0836 Fieldin(end,:)=Fieldin(2,:);
0837
0838 disp(['Interpolating U:',num2str(k),'...'])
0839 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0840 Fieldout=inpaint_nans(Field0,0);
0841 disp('Done.')
0842
0843 disp(['Zeroing U:',num2str(k),' within topography'])
0844 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0845 disp('Done.')
0846 end
36eeb8dc0c Bayl*0847
b9c44f1d6b Bayl*0848 for j=1:length(pickout)
0849 fout=netcdf([dirout '/' pickout(j).name],'write');
0850 Xhere=fout{'Xp1'}(:);
0851 Yhere=fout{'Y'}(:);
0852 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0853 fout{'U'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0854 close(fout)
0855 end
0856
36eeb8dc0c Bayl*0857 %%%%%%%%%%%%%% gUnm1
0858
0859 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0860 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0861 Fieldin=NaN*ones(size(xcin));
0862 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0863
b9c44f1d6b Bayl*0864 for j=1:length(pickin)
0865 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0866 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0867 Xhere=fin{'Xp1'}(:);
0868 Yhere=fin{'Y'}(:);
36eeb8dc0c Bayl*0869 Fieldinhere=squeeze(fin{'gUnm1'}(snap,k,:,:));
0870 HFacinhere=squeeze(gin{'HFacW'}(k,:,:));
b9c44f1d6b Bayl*0871 for ii=1:length(Xhere)
0872 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0873 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0874 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0875 if (HFacinhere(jj,ii)==0)
0876 Fieldin(jjj,iii)=NaN;
0877 end
b9c44f1d6b Bayl*0878 end
0879 end
0880 close(fin)
36eeb8dc0c Bayl*0881 close(gin)
b9c44f1d6b Bayl*0882 end
36eeb8dc0c Bayl*0883 % This utilizes periodic geometry
0884 %Note that there is already a periodic repeat in gUnm1...
b9c44f1d6b Bayl*0885 Fieldin(:,1)=Fieldin(:,end-2);
0886 Fieldin(:,end)=Fieldin(:,3);
0887
36eeb8dc0c Bayl*0888 Fieldin(1,:)=Fieldin(end-1,:);
0889 Fieldin(end,:)=Fieldin(2,:);
0890
0891 disp(['Interpolating gUnm1:',num2str(k),'...'])
0892 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0893 Fieldout=inpaint_nans(Field0,0);
0894 disp('Done.')
0895
0896 disp(['Zeroing gUnm1:',num2str(k),' within topography'])
0897 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0898 disp('Done.')
0899 end
36eeb8dc0c Bayl*0900
b9c44f1d6b Bayl*0901 for j=1:length(pickout)
0902 fout=netcdf([dirout '/' pickout(j).name],'write');
0903 Xhere=fout{'Xp1'}(:);
0904 Yhere=fout{'Y'}(:);
0905 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0906 fout{'gUnm1'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0907 close(fout)
0908 end
0909
0910 %V,gVnm1 is on XV,Rc
0911
36eeb8dc0c Bayl*0912 %Note that there is already a periodic repeat in Vin...
0913 Yp1inext=[2*Yp1in(1)-Yp1in(2);Yp1in;2*Yp1in(end)-Yp1in(end-1)];
0914
0915 [xcin,ycin]=meshgrid(Xinext,Yp1inext);
b9c44f1d6b Bayl*0916 [xcout,ycout]=meshgrid(Xout,Yp1out);
0917
36eeb8dc0c Bayl*0918 %%%%%%%%%%%%% HFacSoutk
0919
0920 HFacout=ones(size(xcout));
0921 HFacoutk=ones([length(Zcomp) size(xcout)]);
0922
0923 disp(['Calculating HFacS...'])
0924 for j=1:length(gridout)
0925 gout=netcdf([dirout '/' gridout(j).name],'nowrite');
0926 Xhere=gout{'X'}(:);
0927 Yhere=gout{'Yp1'}(:);
0928 HFacouthere=gout{'HFacS'}(:,:,:);
0929 for ii=1:length(Xhere)
0930 for jj=1:length(Yhere)
0931 [jjj,iii]=find((Xhere(ii)==xcout).*(Yhere(jj)==ycout));
0932 HFacoutk(:,jjj,iii)=HFacouthere(:,jj,ii);
0933 end
0934 end
0935 close(gout)
0936 end
0937 clear HFacouthere
0938 disp(['Done.'])
0939
0940 %%%%%%%%%%%%%% V
0941
0942 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0943 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0944 Fieldin=NaN*ones(size(xcin));
0945 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0946
b9c44f1d6b Bayl*0947 for j=1:length(pickin)
0948 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*0949 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*0950 Xhere=fin{'X'}(:);
0951 Yhere=fin{'Yp1'}(:);
36eeb8dc0c Bayl*0952 Fieldinhere=squeeze(fin{'V'}(snap,k,:,:));
0953 HFacinhere=squeeze(gin{'HFacS'}(k,:,:));
b9c44f1d6b Bayl*0954 for ii=1:length(Xhere)
0955 for jj=1:length(Yhere)
36eeb8dc0c Bayl*0956 [jjj,iii]=find((Xhere(ii)==xcin).*(Yhere(jj)==ycin));
0957 Fieldin(jjj,iii)=Fieldinhere(jj,ii);
0958 if (HFacinhere(jj,ii)==0)
0959 Fieldin(jjj,iii)=NaN;
0960 end
b9c44f1d6b Bayl*0961 end
0962 end
0963 close(fin)
36eeb8dc0c Bayl*0964 close(gin)
b9c44f1d6b Bayl*0965 end
36eeb8dc0c Bayl*0966 % This utilizes periodic geometry
0967 %Note that there is already a periodic repeat in V...
0968 Fieldin(:,1)=Fieldin(:,end-1);
0969 Fieldin(:,end)=Fieldin(:,2);
b9c44f1d6b Bayl*0970
36eeb8dc0c Bayl*0971 Fieldin(1,:)=Fieldin(end-2,:);
0972 Fieldin(end,:)=Fieldin(3,:);
0973
0974 disp(['Interpolating V:',num2str(k),'...'])
0975 Field0=griddata(xcin,ycin,Fieldin,xcout,ycout,'linear',{'Qt','Qbb','Qc','Qz'});
0976 Fieldout=inpaint_nans(Field0,0);
0977 disp('Done.')
0978
0979 disp(['Zeroing V:',num2str(k),' within topography'])
0980 Fieldoutk(k,:,:)=Fieldout.*squeeze(HFacoutk(k,:,:));
b9c44f1d6b Bayl*0981 disp('Done.')
0982 end
36eeb8dc0c Bayl*0983
b9c44f1d6b Bayl*0984 for j=1:length(pickout)
0985 fout=netcdf([dirout '/' pickout(j).name],'write');
0986 Xhere=fout{'X'}(:);
0987 Yhere=fout{'Yp1'}(:);
0988 [jj,ii]=find((Xhere(1)<=xcout).*(xcout<=Xhere(end)).*(Yhere(1)<=ycout).*(ycout<=Yhere(end)));
0989 fout{'V'}(:,:,:)=Fieldoutk(:,min(jj):max(jj),min(ii):max(ii));
0990 close(fout)
0991 end
0992
36eeb8dc0c Bayl*0993 %%%%%%%%%%%%%% gVnm1
0994
0995 Fieldoutk=zeros([length(Zcomp) size(xcout)]);
0996 for k=1:length(Zcomp)
b9c44f1d6b Bayl*0997 Fieldin=NaN*ones(size(xcin));
0998 Fieldout=NaN*ones(size(xcout));
36eeb8dc0c Bayl*0999
b9c44f1d6b Bayl*1000 for j=1:length(pickin)
1001 fin=netcdf([dirin '/' pickin(j).name],'nowrite');
36eeb8dc0c Bayl*1002 gin=netcdf([dirin '/' gridin(j).name],'nowrite');
b9c44f1d6b Bayl*1003 Xhere=fin{'X'}(:);
1004 Yhere=fin{'Yp1'}(:);
36eeb8dc0c Bayl*1005 Fieldinhere=squeeze(fin{'gVnm1'}(snap,k,:,:));
1006 HFacinhere=squeeze(gin{'HFacS'}(k,:,:));
b9c44f1d6b Bayl*1007 for ii=1:length(Xhere)
1008 for jj=1:length(Yhere)
36eeb8dc0c Bayl*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
b9c44f1d6b Bayl*1014 end
1015 end
1016 close(fin)
36eeb8dc0c Bayl*1017 close(gin)
b9c44f1d6b Bayl*1018 end
36eeb8dc0c Bayl*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);
b9c44f1d6b Bayl*1023
36eeb8dc0c Bayl*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,:,:));
b9c44f1d6b Bayl*1034 disp('Done.')
1035 end
36eeb8dc0c Bayl*1036
b9c44f1d6b Bayl*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
