Warning, /verification/fizhi-cs-aqualev20/scripts/assemble_ML.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
7ae8fb32b5 Andr*0001 %
0002 %
0003 % Ed Hill
0004 %
0005 % Generate the APE ML01--ML07 and ME01--ME07 fields:
0006 %
0007 % For all time steps:
0008 % For all variables:
0009 % For all facets:
0010 % interpolate to the P17 pressure levels:
0011 % compute the "derived" fields:
0012 % sum:
0013 % Compute the averages from the sums:
0014 % Rotate vector components as necessary:
0015 % Convert units as necessary:
0016 % Regrid onto lat-lon:
0017 % Write netCDF output with all attributes:
0018 %
0019
0020 %======================================================================
0021 %
0022 % Define the connections between the APE "ML" variables and the
0023 % MITgcm diagnostics output
0024 %
0025 nfacets = 6;
0026 % bname = 'data/MLave.t%03d.nc';
0027 bname = 'data/ML_ave_%d.nc';
0028 gname = 'data/grid.t%03d.nc';
0029 r_con_name = 'scrip_regrid/rmp_CS32_to_LL128x64_conserv.nc';
0030 r_dwt_name = 'scrip_regrid/rmp_CS32_to_LL128x64_distwgt.nc';
0031 oname = { 'ML_fields.nc' };
0032
0033 flags.regrid_type_for_vec = 2;
0034 flags.debug_lev = 0;
0035 flags.punits_conv = 0.01;
0036
0037 % Get the following from:
0038 % ! (cd ../../ape_data_specs/ ; ./ml_parse.sh)
0039 vars = {};
0040 vars{1} = {'ML01','ml_u','UVEL','eastward_wind','m s-1','u','1','ee'};
0041 vars{2} = {'ML02','ml_v','VVEL','northward_wind','m s-1','v','1','ee'};
0042 vars{3} = {'ML03','ml_t','ML03','air_temperature','K','-','1','ec'};
0043 vars{4} = {'ML04','ml_om','WVEL','omega','Pa s-1','-','1','ee'};
0044 vars{5} = {'ML05','ml_phi','ML05','geopotential_height','m','-','1','ee'};
0045 vars{6} = {'ML06','ml_q','SALT','specific_humidity','kg kg-1','-','1','cc'};
0046 vars{7} = {'ML07','ml_rh','RELHUM','relative_humidity','percent','-','1','cc'};
0047 vars{8} = {'ML08','ml_th','THETA','potential_temperature','K','-','1','ec'};
0048
0049 % Geopotential reference values
0050 % load geopot_20 ; geopot = geopot_20(:,2); clear geopot_20
0051 phiref = ...
0052 [ -0.739878953 646.302002 1338.38696 2894.67627 4099.63135 ...
0053 5484.93359 7116.62549 9115.08008 10321.4688 11741.3574 13494.4102 ...
0054 15862.4404 17833.5605 19667.8887 22136.1934 23854.2266 26366.9375 ];
0055
0056 %======================================================================
0057 %
0058 % Open the input files and get the (minimum) number of time steps
0059 % available across all the files
0060 %
0061 disp('Finding available time steps')
0062
0063 ncl = {};
0064 nTmin = 0;
0065 ifacet = 1;
0066 for ifacet = 1:nfacets
0067 fname = sprintf(bname,ifacet);
0068 nc = netcdf(fname, 'nowrite');
0069 ncl{ifacet} = nc;
0070 nT = length( nc('T') );
0071 if (ifacet == 1)
0072 nTmin = nT;
0073 else
0074 nTmin = min(nTmin, nT);
0075 end
0076 nc = [];
0077 end
0078 disp(sprintf(' total iterations found = %d',nTmin));
0079
0080
0081 %======================================================================
0082 %
0083 % Interpolate to the reference pressure levels and then compute the
0084 % derived fields and the temporal sums of all fields on the original
0085 % MITgcm grid
0086 %
0087 disp('Computing Sums')
0088
0089 Pref = [ 1 2 3 5 7 10 15 20 25 30 40 50 60 70 85 92.5 100 ] * 10;
0090 Pref = Pref(length(Pref):-1:1);
0091
0092 dat = [];
0093 it = 1;
0094 iv = 1;
0095 ifacet = 1;
0096 it0 = 1;
0097 for it = [ 1 ] % months.i_start:min(nTmin,months.i_end)
0098 disp(sprintf(' iteration = %d',it));
0099
0100 for iv = 1:length(vars)
0101 for ifacet = 1:nfacets
0102
0103 t0 = [];
0104 t1 = [];
0105 vtmp = [];
0106
0107 % cull unknown variables
0108 if ( vars{iv}{3}(1) == '-' )
0109 continue
0110 end
0111
0112 % skip over computed quantities
0113 if strcmp( vars{iv}{1} , vars{iv}{3} )
0114 continue
0115 end
0116
0117 % WARN about missing variables
0118 nc = ncl{ifacet};
0119 if prod(size(nc{vars{iv}{3}})) == 0
0120 if it == 1
0121 str = 'var "%s" does not exist in file "%s"';
0122 disp(sprintf([' Warning: ' str], vars{iv}{3}, name(nc)));
0123 end
0124 continue
0125 end
0126
0127 % get the data
0128 t0 = squeeze( nc{ vars{iv}{3} }(it,:,:,:) );
0129
0130 % horizontally interpolate to u,v to mass points
0131 %
0132 % this has changed so that now the "+1" points are garbage
0133 % and must be ignored
0134 switch vars{iv}{6}
0135 case 'u'
0136 t1 = t0(:,:,1:32);
0137 case 'v'
0138 t1 = t0(:,1:32,:);
0139 otherwise
0140 t1 = t0;
0141 end
0142
0143 vtmp = t1;
0144
0145 % sum for averages
0146 if it0 == 1
0147 dat(iv,ifacet).n = 1;
0148 dat(iv,ifacet).a = vtmp;
0149 else
0150 dat(iv,ifacet).n = dat(iv,ifacet).n + 1;
0151 dat(iv,ifacet).a = dat(iv,ifacet).a + vtmp;
0152 end
0153
0154 end
0155 end
0156
0157 % fill in the derived fields
0158 for iv = 1:length(vars)
0159 for ifacet = 1:nfacets
0160
0161 tmp = [];
0162 vtmp = [];
0163 nc = ncl{ifacet};
0164
0165 if strcmp(vars{iv}{1},'ML03')
0166 % ML03 = THETA * (P/P0)^0.286
0167 tmp = squeeze( nc{ 'THETA' }(it,:,:,:) );
0168 for ii = 1:32
0169 for jj = 1:32
0170 tmp(:,ii,jj) = tmp(:,ii,jj) .* ( Pref' ./ 1000 ).^0.286;
0171 end
0172 end
0173 end
0174
0175 if strcmp(vars{iv}{1},'ML05')
0176 % ML05 = (PHIHYD/9.81) + ref
0177 tmp = squeeze( nc{ 'PHIHYD' }(it,:,:,:) );
0178 nn = size(tmp);
0179 for ii = 1:32
0180 for jj = 1:32
0181 tmp(:,ii,jj) = squeeze(tmp(:,ii,jj))/9.81 + phiref'; % + geopot;
0182 end
0183 end
0184 end
0185
0186 if prod(size(tmp)) > 0
0187
0188 if it0 == 1
0189 dat(iv,ifacet).n = 1;
0190 dat(iv,ifacet).a = tmp;
0191 else
0192 dat(iv,ifacet).n = dat(iv,ifacet).n + 1;
0193 dat(iv,ifacet).a = dat(iv,ifacet).a + tmp;
0194 end
0195 end
0196
0197 end
0198 end
0199
0200 it0 = 0;
0201
0202 end
0203 clear tmp
0204
0205 % close the input files
0206 for ifacet = 1:nfacets
0207 nc = close( ncl{ifacet} );
0208 end
0209 clear ncl;
0210
0211 if 1 == 0
0212 for iv = 1:length(vars)
0213 for ifacet = 1:nfacets
0214 for iz = 1:5:17
0215 aa = squeeze( dat(iv,ifacet).a(iz,:,:) );
0216 surf(aa), view(2), shading flat, colorbar
0217 disp([iv ifacet iz]);
0218 pause(1)
0219 end
0220 end
0221 end
0222 end
0223
0224 % interpolate u,v to mass points
0225 nc = 32;
0226 ncp = nc + 1;
0227
0228 for iv = 1:length(vars)
0229 switch vars{iv}{6}
0230 case 'u'
0231
0232 nr = size(dat(iv,1).a,1);
0233 u3d = zeros(nc,nc,nr,6);
0234 v3d = zeros(nc,nc,nr,6);
0235 for fi = 1:6
0236 u3d(:,:,:,fi) = permute( dat(iv, fi).a, [3 2 1] );
0237 v3d(:,:,:,fi) = permute( dat(iv+1,fi).a, [3 2 1] );
0238 end
0239
0240 u6t=zeros(ncp,nc,nr,6);
0241 v6t=zeros(nc,ncp,nr,6);
0242 u6t([1:nc],:,:,:)=u3d(:,:,:,:);
0243 v6t(:,[1:nc],:,:)=v3d(:,:,:,:);
0244
0245 u6t(ncp,[1:nc],:,1)=u3d(1,[1:nc],:,2);
0246 u6t(ncp,[1:nc],:,2)=v3d([nc:-1:1],1,:,4);
0247 u6t(ncp,[1:nc],:,3)=u3d(1,[1:nc],:,4);
0248 u6t(ncp,[1:nc],:,4)=v3d([nc:-1:1],1,:,6);
0249 u6t(ncp,[1:nc],:,5)=u3d(1,[1:nc],:,6);
0250 u6t(ncp,[1:nc],:,6)=v3d([nc:-1:1],1,:,2);
0251
0252 v6t([1:nc],ncp,:,1)=u3d(1,[nc:-1:1],:,3);
0253 v6t([1:nc],ncp,:,2)=v3d([1:nc],1,:,3);
0254 v6t([1:nc],ncp,:,3)=u3d(1,[nc:-1:1],:,5);
0255 v6t([1:nc],ncp,:,4)=v3d([1:nc],1,:,5);
0256 v6t([1:nc],ncp,:,5)=u3d(1,[nc:-1:1],:,1);
0257 v6t([1:nc],ncp,:,6)=v3d([1:nc],1,:,1);
0258
0259 for fi = 1:6
0260 dat(iv, fi).a = permute( ...
0261 0.5*( u6t([1:nc],:,:,fi) + u6t([2:ncp],:,:,fi) ), [3 2 1] );
0262 dat(iv+1,fi).a = permute( ...
0263 0.5*( v6t(:,[1:nc],:,fi) + v6t(:,[2:ncp],:,fi) ), [3 2 1] );
0264 end
0265
0266 end
0267 end
0268
0269 %======================================================================
0270 %
0271 % Compute the averages from the temporal sums
0272 %
0273 disp('Computing Averages from Sums')
0274
0275 ave = [];
0276 for iv = 1:length(vars)
0277 ind = length(ave) + 1;
0278 for ifacet = 1:nfacets
0279 if length(dat(iv,ifacet).n) == 0
0280 continue
0281 end
0282 ave(ind).v(:,:,:,ifacet) = dat(iv,ifacet).a ./ dat(iv,ifacet).n;
0283 ave(ind).ivar = iv;
0284 end
0285
0286 if 1 == 0
0287 aa = squeeze(ave(ind).v(1,:,:,1));
0288 surf(aa), view(2), shading interp, colorbar, pause(3)
0289 close all
0290 end
0291
0292 end
0293
0294 clear dat;
0295
0296 %======================================================================
0297 %
0298 % Rotate vector components & convert units
0299 %
0300 clear ig;
0301 gvars = { 'XC','YC','dxF','dyF','rA','XG','YG','dxV', ...
0302 'dyU','rAz','dxC','dyC','rAw','rAs','dxG','dyG' };
0303 for iface = 1:6
0304 fname = sprintf(gname, iface);
0305 nc = netcdf(fname, 'nowrite');
0306 if iface == 1
0307 ig.ne = length(nc('X'));
0308 end
0309 for jj = 1:length(gvars)
0310 comm = [ 'ig.' gvars{jj} '(:,:,iface) = nc{''' ...
0311 gvars{jj} '''}(:);'];
0312 eval(comm);
0313 end
0314 nc = close(nc);
0315 end
0316
0317 ne = ig.ne;
0318 n1 = ne - 1;
0319 dux = zeros(size(ig.XC));
0320 duy = zeros(size(ig.XC));
0321 dvx = zeros(size(ig.XC));
0322 dvy = zeros(size(ig.XC));
0323 % dux(:,:,:) = diff(ig.XG(:,1:ne,:),1,1);
0324 % dvx(:,:,:) = diff(ig.XG(1:ne,:,:),1,2);
0325 % duy(:,:,:) = diff(ig.YG(:,1:ne,:),1,1);
0326 % dvy(:,:,:) = diff(ig.YG(1:ne,:,:),1,2);
0327 %
0328 % Note the first two dimensions below are permuted relative to the
0329 % ordering in the *.mitgrid files (commented-out above).
0330 dux(:,:,:) = diff(ig.XG(1:ne,:,:),1,2);
0331 dvx(:,:,:) = diff(ig.XG(:,1:ne,:),1,1);
0332 duy(:,:,:) = diff(ig.YG(1:ne,:,:),1,2);
0333 dvy(:,:,:) = diff(ig.YG(:,1:ne,:),1,1);
0334 dux = dux + 360*double(dux < 180);
0335 dux = dux - 360*double(dux > 180); % squeeze([ min(min(dux)) max(max(dux)) ])
0336 duy = duy + 360*double(duy < 180);
0337 duy = duy - 360*double(duy > 180); % squeeze([ min(min(duy)) max(max(duy)) ])
0338 dvx = dvx + 360*double(dvx < 180);
0339 dvx = dvx - 360*double(dvx > 180); % squeeze([ min(min(dvx)) max(max(dvx)) ])
0340 dvy = dvy + 360*double(dvy < 180);
0341 dvy = dvy - 360*double(dvy > 180); % squeeze([ min(min(dvy)) max(max(dvy)) ])
0342 dut = sqrt(dux.^2 + duy.^2); % squeeze([ min(min(dut)) max(max(dut)) ])
0343 dvt = sqrt(dvx.^2 + dvy.^2); % squeeze([ min(min(dvt)) max(max(dvt)) ])
0344 ig.llux = dux ./ dut; % squeeze([ min(min(ig.llux)) max(max(ig.llux)) ])
0345 ig.lluy = duy ./ dut; % squeeze([ min(min(ig.lluy)) max(max(ig.lluy)) ])
0346 ig.llvx = dvx ./ dvt; % squeeze([ min(min(ig.llvx)) max(max(ig.llvx)) ])
0347 ig.llvy = dvy ./ dvt; % squeeze([ min(min(ig.llvy)) max(max(ig.llvy)) ])
0348 clear dux duy dvx dvy dut dvt ne n1
0349
0350 if flags.regrid_type_for_vec ~= 2
0351 % Vector rotation
0352 for ia = 1:length(ave)
0353 if strcmp(vars{ ave(ia).ivar }{6}, 'u')
0354 % llu = u .* llux + v .* llvx;
0355 % llv = u .* lluy + v .* llvy;
0356 for iz = 1:size(ave(ia).v,1)
0357 llu = squeeze( ave(ia).v(iz,:,:,:) ) .* ig.llux ...
0358 + squeeze( ave(ia+1).v(iz,:,:,:) ) .* ig.llvx;
0359 llv = squeeze( ave(ia).v(iz,:,:,:) ) .* ig.lluy ...
0360 + squeeze( ave(ia+1).v(iz,:,:,:) ) .* ig.llvy;
0361 ave(ia).v(iz,:,:,:) = llu;
0362 ave(ia+1).v(iz,:,:,:) = llv;
0363 end
0364 end
0365 end
0366 end
0367 clear llu llv
0368
0369 % Units conversion
0370 for ia = 1:length(ave)
0371 if strcmp(vars{ ave(ia).ivar }{7}, '-') ~= 1
0372 eval(['fac = ' vars{ ave(ia).ivar }{7} ';']);
0373 if fac ~= 1
0374 ave(ia).v = fac .* ave(ia).v;
0375 end
0376 else
0377
0378 end
0379 end
0380
0381 %======================================================================
0382 %
0383 % Regrid
0384 %
0385 % JMC suggested an evenly spaced Lat-Lon at: 128x64
0386 %
0387 disp('Regridding')
0388
0389 og = [];
0390 og.nlat = 64;
0391 og.nlon = 128;
0392 og.latcell = 180/og.nlat;
0393 og.loncell = 360/og.nlon;
0394 og.lat = linspace(-90+og.latcell/2, 90-og.latcell/2, og.nlat);
0395 og.lon = linspace(0+og.loncell/2, 360-og.loncell/2, og.nlon);
0396 og.latbnd(:,1) = og.lat - og.latcell/2.0;
0397 og.latbnd(:,2) = og.lat + og.latcell/2.0;
0398 og.lonbnd(:,1) = og.lon - og.loncell/2.0;
0399 og.lonbnd(:,2) = og.lon + og.loncell/2.0;
0400
0401
0402 rgvar = { 'src_address', 'dst_address', 'remap_matrix' };
0403 rtype = { { 'con', r_con_name}, {'dwt', r_dwt_name } };
0404 for kk = 1:length(rtype)
0405 nc = netcdf(rtype{kk}{2}, 'nowrite');
0406 for jj = 1:length(rgvar)
0407 comm = sprintf('rgrid.%s.%s = nc{''%s''}(:);',...
0408 rtype{kk}{1},rgvar{jj},rgvar{jj});
0409 eval(comm);
0410 % disp(comm);
0411 end
0412 nc = close(nc);
0413 end
0414
0415 if ( flags.regrid_type_for_vec > 0 )
0416 % >> x=rdmds('XC');
0417 % >> y=rdmds('YC');
0418 % >> t=rdmds('Ttave.0000513360');
0419 % >> xi=-179:2:180;yi=-89:2:90;
0420 % >> del=cube2latlon_preprocess(x,y,xi,yi);
0421 % >> ti=cube2latlon_fast(del,t);
0422 aa_regrid.x = rdmds('aa_regrid/XC');
0423 aa_regrid.y = rdmds('aa_regrid/YC');
0424 % u = rdmds('aa_regrid/U.0000000000');
0425 % v = rdmds('aa_regrid/V.0000000000');
0426 % del = cube2latlon_preprocess(LON,LAT,xc,yc);
0427 %
0428 % aa_regrid.del = ...
0429 % cube2latlon_preprocess(aa_regrid.x,aa_regrid.y, ...
0430 % og.lon,og.lat);
0431 end
0432
0433 % save ML_ave
0434
0435 ii = 1;
0436 for ii = 1:length(ave)
0437
0438 for iz = 1:size( ave(ii).v, 1 )
0439
0440 siz = size( ave(ii).v );
0441 dcp = zeros(prod(siz(2:4)),1);
0442 nn = 0;
0443 for fk = 1:size( ave(ii).v, 4 )
0444 for jk = 1:size( ave(ii).v, 3 )
0445 for ik = 1:size( ave(ii).v, 2 )
0446 nn = nn + 1;
0447 dcp(nn) = ave(ii).v(iz,ik,jk,fk);
0448 end
0449 end
0450 end
0451 lld = zeros(og.nlat*og.nlon,1);
0452
0453 if ( ( flags.regrid_type_for_vec > 0 ) ...
0454 && ( vars{ave(ii).ivar}{6} == 'u' ...
0455 || vars{ave(ii).ivar}{6} == 'v' ) )
0456
0457 if flags.regrid_type_for_vec == 1
0458
0459 % use distributed weights for vectors
0460 for jj = 1:length(rgrid.dwt.src_address)
0461 lld(rgrid.dwt.dst_address(jj)) = ...
0462 lld(rgrid.dwt.dst_address(jj)) ...
0463 + ( dcp(rgrid.dwt.src_address(jj)) ...
0464 * rgrid.dwt.remap_matrix(jj,1) );
0465 end
0466 ave(ii).r(iz,:,:) = reshape(lld',og.nlat,og.nlon);
0467
0468 elseif ( flags.regrid_type_for_vec == 2 ...
0469 && vars{ave(ii).ivar}{6} == 'u' )
0470
0471 % use Alistair's uvcube2latlon() for vectors
0472 nsiz = size(ave(ii).v);
0473 by6u = zeros([ nsiz(2)*nsiz(4) nsiz(3) ]);
0474 by6v = zeros([ nsiz(2)*nsiz(4) nsiz(3) ]);
0475 for kk = 1:size( ave(ii).v, 4 )
0476 is = 1 + (kk-1)*nsiz(2);
0477 ie = is + nsiz(2) - 1;
0478 by6u(is:ie,:) = squeeze( ave(ii ).v(iz,:,:,kk) )';
0479 by6v(is:ie,:) = squeeze( ave(ii+1).v(iz,:,:,kk) )';
0480 end
0481 % xi=-179:2:180; yi=-89:2:90;
0482 xi = og.lon;
0483 ind = find(xi > 179.5);
0484 xi(ind) = xi(ind) - 360;
0485 yi = og.lat;
0486 [ul,vl] = uvacube2latlon(aa_regrid.x,aa_regrid.y, ...
0487 by6u,by6v, xi,yi);
0488 ave(ii ).r(iz,:,:) = ul';
0489 ave(ii+1).r(iz,:,:) = vl';
0490
0491 if flags.debug_lev > 0
0492 figure(1), subplot(1,1,1)
0493 subplot(2,1,1), surf(ul'), axis equal, view(2), shading flat
0494 subplot(2,1,2), surf(vl'), axis equal, view(2), shading flat
0495 end
0496
0497 elseif ( vars{ave(ii).ivar}{6} == 'v' )
0498
0499 if iz == 1
0500 disp(sprintf( ' %s was computed along with %s', ...
0501 vars{ave(ii).ivar}{2}, vars{ave(ii-1).ivar}{2} ));
0502 end
0503
0504 elseif ( vars{ave(ii).ivar}{6} == 'u' )
0505
0506 disp('Error: please specify regrid type for vectors');
0507 disp(' using "flags.regrid_type_for_vec"');
0508
0509 end
0510 else
0511 % conservative for scalars
0512 for jj = 1:length(rgrid.con.src_address)
0513 lld(rgrid.con.dst_address(jj)) = ...
0514 lld(rgrid.con.dst_address(jj)) ...
0515 + ( dcp(rgrid.con.src_address(jj)) ...
0516 * rgrid.con.remap_matrix(jj,1) );
0517 end
0518 ave(ii).r(iz,:,:) = reshape(lld',og.nlat,og.nlon);
0519 end
0520
0521 end
0522
0523 end
0524
0525
0526 %======================================================================
0527 %
0528 % Write netCDF output
0529 %
0530 disp('Writing netCDF output')
0531
0532 nc = netcdf(oname{1}, 'clobber');
0533 nc.title = [ 'Aqua Planet: ' ...
0534 'Single-Level 2-D Means from "Example" Experiment' ];
0535 nc.institution = 'MIT Dept. of EAPS, Cambridge, MA, USA';
0536 nc.source = [ 'MITgcm: ' ];
0537 nc.Conventions = 'CF-1.0';
0538 nc.history = [ 'Original data produced: ' '2002/08/20' ];
0539
0540 nc('pres') = length(Pref);
0541 nc('lon') = og.nlon;
0542 nc('lat') = og.nlat;
0543 nc('bnd') = 2;
0544
0545 nc{'pres'} = ncdouble('pres');
0546 nc{'pres'}.standard_name = 'air_pressure';
0547 nc{'pres'}.units = 'Pa';
0548 nc{'pres'}(:) = Pref*100;
0549
0550 nc{'lon'} = ncdouble('lon');
0551 nc{'lon'}.standard_name = 'longitude';
0552 nc{'lon'}.units = 'degrees_east';
0553 nc{'lon'}.bounds = 'lon_bnds';
0554 nc{'lon'}(:) = og.lon;
0555
0556 nc{'lat'} = ncdouble('lat');
0557 nc{'lat'}.standard_name = 'latitude';
0558 nc{'lat'}.units = 'degrees_north';
0559 nc{'lat'}.bounds = 'lat_bnds';
0560 nc{'lat'}(:) = og.lat;
0561
0562 nc{'lon_bnds'} = ncdouble('lon','bnd');
0563 nc{'lon_bnds'}.ape_name = 'longitude cell bounds';
0564 nc{'lon_bnds'}.units = 'degrees_east';
0565 nc{'lon_bnds'}(:) = og.lonbnd;
0566
0567 nc{'lat_bnds'} = ncdouble('lat','bnd');
0568 nc{'lat_bnds'}.ape_name = 'latitude cell bounds';
0569 nc{'lat_bnds'}.units = 'degrees_north';
0570 nc{'lat_bnds'}(:) = og.latbnd;
0571
0572 for ii = 1:length(ave)
0573
0574 iv = ave(ii).ivar;
0575 nc{ vars{iv}{2} } = ncfloat( 'pres', 'lat', 'lon' );
0576 nc{ vars{iv}{2} }.ape_name = vars{iv}{4};
0577 nc{ vars{iv}{2} }.units = vars{iv}{5};
0578 nc{ vars{iv}{2} }.FillValue_ = 1.0e20;
0579
0580 nc{ vars{iv}{2} }(:) = ave(ii).r;
0581
0582 end
0583
0584 nc = close(nc);
