Warning, /verification/fizhi-cs-aqualev20/scripts/assemble_TR.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 TR fields:
0006 %
0007
0008 %======================================================================
0009 %
0010 % Define the connections between the APE "ML" variables and the
0011 % MITgcm diagnostics output
0012 %
0013 nfacets = 6;
0014 months.i_start = 7;
0015 months.i_end = 42;
0016 bname = 'data/TRall.t%03d.nc';
0017 gname = 'data/grid.t%03d.nc';
0018 r_con_name = 'scrip_regrid/rmp_CS32_to_LL128x64_conserv.nc';
0019 r_dwt_name = 'scrip_regrid/rmp_CS32_to_LL128x64_distwgt.nc';
0020 oname = { 'TR_fields.nc' };
0021
0022 flags.regrid_type_for_vec = 2;
0023 flags.debug_lev = 0;
0024 flags.punits_conv = 0.01;
0025
0026 % Get the following from:
0027 % ! (cd ../../ape_data_specs/ ; ./tr_parse.sh)
0028 vars = {};
0029 vars{1} = {'TR01','tr_tppn','PREACC','precipitation_flux','kg m-2 s-1','-','0.000011574074','ee'};
0030 vars{2} = {'TR02','tr_lw_toa','OLR','toa_outgoing_longwave_flux','W m-2','-','1','cc'};
0031 vars{3} = {'TR03','tr_om500','WVEL','omega','Pa s-1','-','1','ec'};
0032 vars{4} = {'TR04','tr_u250','UVEL','eastward_wind','m s-1','u','1','ee'};
0033 vars{5} = {'TR05','tr_v250','VVEL','northward_wind','m s-1','v','1','ee'};
0034 vars{6} = {'TR06','tr_mslp','RSURF','air_pressure_at_sea_level','Pa','-','1','cc'};
0035
0036
0037 %======================================================================
0038 %
0039 % Open the input files and get the (minimum) number of time steps
0040 % available across all the files
0041 %
0042 disp('Finding available time steps')
0043
0044 ncl = {};
0045 nTmin = 0;
0046 ifacet = 1;
0047 for ifacet = 1:nfacets
0048 fname = sprintf(bname,ifacet);
0049 nc = netcdf(fname, 'nowrite');
0050 ncl{ifacet} = nc;
0051 T = nc{'T'}(:);
0052 nT = length(T);
0053 if (ifacet == 1)
0054 nTmin = nT;
0055 else
0056 nTmin = min(nTmin, nT);
0057 end
0058 nc = [];
0059 end
0060 disp(sprintf(' total iterations found = %d',nTmin));
0061
0062
0063 %======================================================================
0064 %
0065 og = [];
0066 og.nlat = 64;
0067 og.nlon = 128;
0068 og.latcell = 180/og.nlat;
0069 og.loncell = 360/og.nlon;
0070 og.lat = linspace(-90+og.latcell/2, 90-og.latcell/2, og.nlat);
0071 og.lon = linspace(0+og.loncell/2, 360-og.loncell/2, og.nlon);
0072 og.latbnd(:,1) = og.lat - og.latcell/2.0;
0073 og.latbnd(:,2) = og.lat + og.latcell/2.0;
0074 og.lonbnd(:,1) = og.lon - og.loncell/2.0;
0075 og.lonbnd(:,2) = og.lon + og.loncell/2.0;
0076
0077 rgvar = { 'src_address', 'dst_address', 'remap_matrix' };
0078 rtype = { { 'con', r_con_name}, {'dwt', r_dwt_name } };
0079 for kk = 1:length(rtype)
0080 nc = netcdf(rtype{kk}{2}, 'nowrite');
0081 for jj = 1:length(rgvar)
0082 comm = sprintf('rgrid.%s.%s = nc{''%s''}(:);',...
0083 rtype{kk}{1},rgvar{jj},rgvar{jj});
0084 eval(comm);
0085 % disp(comm);
0086 end
0087 nc = close(nc);
0088 end
0089
0090 aa_regrid.x = rdmds('aa_regrid/XC');
0091 aa_regrid.y = rdmds('aa_regrid/YC');
0092
0093 disp('Reading -- Regridding -- Writing')
0094
0095 nc_outf = netcdf(oname{1}, 'clobber');
0096 nc_outf.title = [ 'Aqua Planet: ' ...
0097 'TR Fields' ];
0098 nc_outf.institution = 'MIT Dept. of EAPS, Cambridge, MA, USA';
0099 nc_outf.source = [ 'MITgcm: ' ];
0100 nc_outf.Conventions = 'CF-1.0';
0101 nc_outf.history = [ 'Original data produced: ' '' ];
0102
0103 nc_outf('time') = 0; % record dim
0104 nc_outf('lon') = og.nlon;
0105 nc_outf('lat') = og.nlat;
0106 nc_outf('bnd') = 2;
0107
0108 nc_outf{'time'} = ncdouble('time');
0109 nc_outf{'time'}.standard_name = 'time';
0110 nc_outf{'time'}.units = 'days since 0000-01-01 00:00';
0111 nc_outf{'time'}.bounds = 'time_bnds';
0112 % nc_outf{'time'}(:) = og.lon;
0113
0114 nc_outf{'lon'} = ncdouble('lon');
0115 nc_outf{'lon'}.standard_name = 'longitude';
0116 nc_outf{'lon'}.units = 'degrees_east';
0117 nc_outf{'lon'}.bounds = 'lon_bnds';
0118 nc_outf{'lon'}(:) = og.lon;
0119
0120 nc_outf{'lat'} = ncdouble('lat');
0121 nc_outf{'lat'}.standard_name = 'latitude';
0122 nc_outf{'lat'}.units = 'degrees_north';
0123 nc_outf{'lat'}.bounds = 'lat_bnds';
0124 nc_outf{'lat'}(:) = og.lat;
0125
0126 nc_outf{'time_bnds'} = ncdouble('time','bnd');
0127 nc_outf{'time_bnds'}.ape_name = 'time interval endpoints';
0128 % nc_outf{'time_bnds'}(:) = og.lonbnd;
0129
0130 nc_outf{'lon_bnds'} = ncdouble('lon','bnd');
0131 nc_outf{'lon_bnds'}.ape_name = 'longitude cell bounds';
0132 nc_outf{'lon_bnds'}.units = 'degrees_east';
0133 nc_outf{'lon_bnds'}(:) = og.lonbnd;
0134
0135 nc_outf{'lat_bnds'} = ncdouble('lat','bnd');
0136 nc_outf{'lat_bnds'}.ape_name = 'latitude cell bounds';
0137 nc_outf{'lat_bnds'}.units = 'degrees_north';
0138 nc_outf{'lat_bnds'}(:) = og.latbnd;
0139
0140
0141 for iv = 1:length(vars)
0142 nc_outf{ vars{iv}{2} } = ncfloat( 'time', 'lat', 'lon' );
0143 nc_outf{ vars{iv}{2} }.ape_name = vars{iv}{4};
0144 nc_outf{ vars{iv}{2} }.units = vars{iv}{5};
0145 nc_outf{ vars{iv}{2} }.FillValue_ = 1.0e20;
0146 end
0147
0148 ne = 32;
0149
0150 dat = [];
0151 it = 1;
0152 iv = 1;
0153 ifacet = 1;
0154 it0 = 2880;
0155 for it = [ (it0+1):nTmin ]
0156
0157 if mod(it-1,50) == 0
0158 disp(sprintf([ ' iteration = %6d ' datestr(now) ],it));
0159 end
0160
0161 for iv = 1:length(vars)
0162 for ifacet = 1:nfacets
0163
0164 t0 = [];
0165 vtmp = [];
0166
0167 % WARN about missing variables
0168 nc = ncl{ifacet};
0169 if prod(size(nc{vars{iv}{3}})) == 0
0170 if it == 1
0171 str = 'var "%s" does not exist in file "%s"';
0172 disp(sprintf([' Warning: ' str], vars{iv}{3}, name(nc)));
0173 end
0174 continue
0175 end
0176
0177 % get the data
0178 t0 = nc{ vars{iv}{3} }(it,:,:);
0179
0180 % horizontally interpolate to u,v to mass points
0181 %
0182 % this has changed so that now the "+1" points are garbage
0183 % and must be ignored
0184 dat(iv,ifacet).a = zeros([1 32 32]);
0185 switch vars{iv}{6}
0186 case 'u'
0187 dat(iv,ifacet).a(1,:,:) = t0(:,1:ne);
0188 case 'v'
0189 dat(iv,ifacet).a(1,:,:) = t0(1:ne,:);
0190 otherwise
0191 dat(iv,ifacet).a(1,:,:) = t0;
0192 end
0193
0194 dat(iv,ifacet).n = 1;
0195
0196 end
0197 end
0198
0199 % interpolate u,v to mass points
0200 nc = ne;
0201 ncp = nc + 1;
0202 for iv = 1:length(vars)
0203 switch vars{iv}{6}
0204 case 'u'
0205
0206 nr = size(dat(iv,1).a,1);
0207 u3d = zeros(nc,nc,nr,6);
0208 v3d = zeros(nc,nc,nr,6);
0209 for fi = 1:6
0210 u3d(:,:,:,fi) = permute( dat(iv, fi).a, [3 2 1] );
0211 v3d(:,:,:,fi) = permute( dat(iv+1,fi).a, [3 2 1] );
0212 end
0213
0214 u6t=zeros(ncp,nc,nr,6);
0215 v6t=zeros(nc,ncp,nr,6);
0216 u6t([1:nc],:,:,:)=u3d(:,:,:,:);
0217 v6t(:,[1:nc],:,:)=v3d(:,:,:,:);
0218
0219 u6t(ncp,[1:nc],:,1)=u3d(1,[1:nc],:,2);
0220 u6t(ncp,[1:nc],:,2)=v3d([nc:-1:1],1,:,4);
0221 u6t(ncp,[1:nc],:,3)=u3d(1,[1:nc],:,4);
0222 u6t(ncp,[1:nc],:,4)=v3d([nc:-1:1],1,:,6);
0223 u6t(ncp,[1:nc],:,5)=u3d(1,[1:nc],:,6);
0224 u6t(ncp,[1:nc],:,6)=v3d([nc:-1:1],1,:,2);
0225
0226 v6t([1:nc],ncp,:,1)=u3d(1,[nc:-1:1],:,3);
0227 v6t([1:nc],ncp,:,2)=v3d([1:nc],1,:,3);
0228 v6t([1:nc],ncp,:,3)=u3d(1,[nc:-1:1],:,5);
0229 v6t([1:nc],ncp,:,4)=v3d([1:nc],1,:,5);
0230 v6t([1:nc],ncp,:,5)=u3d(1,[nc:-1:1],:,1);
0231 v6t([1:nc],ncp,:,6)=v3d([1:nc],1,:,1);
0232
0233 for fi = 1:6
0234 dat(iv, fi).a = permute( ...
0235 0.5*( u6t([1:nc],:,:,fi) + u6t([2:ncp],:,:,fi) ), [3 2 1] );
0236 dat(iv+1,fi).a = permute( ...
0237 0.5*( v6t(:,[1:nc],:,fi) + v6t(:,[2:ncp],:,fi) ), [3 2 1] );
0238 end
0239
0240 end
0241 end
0242
0243 % Units conversion
0244 ave = [];
0245 for iv = 1:length(vars)
0246 ind = length(ave) + 1;
0247 for ifacet = 1:nfacets
0248 ave(ind).v(1,:,:,ifacet) = dat(iv,ifacet).a;
0249 ave(ind).ivar = iv;
0250 end
0251 end
0252 for ia = 1:length(ave)
0253 if strcmp(vars{ ave(ia).ivar }{7}, '-') ~= 1
0254 eval(['fac = ' vars{ ave(ia).ivar }{7} ';']);
0255 if fac ~= 1
0256 ave(ia).v = fac .* ave(ia).v;
0257 end
0258 end
0259 end
0260
0261 % regrid
0262 ii = 1;
0263 for ii = 1:length(ave)
0264
0265 for iz = 1:size( ave(ii).v, 1 )
0266
0267 siz = size( ave(ii).v );
0268 dcp = zeros(prod(siz(2:4)),1);
0269 nn = 0;
0270 for fk = 1:size( ave(ii).v, 4 )
0271 for jk = 1:size( ave(ii).v, 3 )
0272 for ik = 1:size( ave(ii).v, 2 )
0273 nn = nn + 1;
0274 dcp(nn) = ave(ii).v(iz,ik,jk,fk);
0275 end
0276 end
0277 end
0278 lld = zeros(og.nlat*og.nlon,1);
0279
0280 if ( ( flags.regrid_type_for_vec > 0 ) ...
0281 && ( vars{ave(ii).ivar}{6} == 'u' ...
0282 || vars{ave(ii).ivar}{6} == 'v' ) )
0283
0284 if flags.regrid_type_for_vec == 1
0285
0286 % use distributed weights for vectors
0287 for jj = 1:length(rgrid.dwt.src_address)
0288 lld(rgrid.dwt.dst_address(jj)) = ...
0289 lld(rgrid.dwt.dst_address(jj)) ...
0290 + ( dcp(rgrid.dwt.src_address(jj)) ...
0291 * rgrid.dwt.remap_matrix(jj,1) );
0292 end
0293 ave(ii).r(iz,:,:) = reshape(lld',og.nlat,og.nlon);
0294
0295 elseif ( flags.regrid_type_for_vec == 2 ...
0296 && vars{ave(ii).ivar}{6} == 'u' )
0297
0298 % use Alistair's uvcube2latlon() for vectors
0299 nsiz = size(ave(ii).v);
0300 by6u = zeros([ nsiz(2)*nsiz(4) nsiz(3) ]);
0301 by6v = zeros([ nsiz(2)*nsiz(4) nsiz(3) ]);
0302 for kk = 1:size( ave(ii).v, 4 )
0303 is = 1 + (kk-1)*nsiz(2);
0304 ie = is + nsiz(2) - 1;
0305 by6u(is:ie,:) = squeeze( ave(ii ).v(iz,:,:,kk) )';
0306 by6v(is:ie,:) = squeeze( ave(ii+1).v(iz,:,:,kk) )';
0307 end
0308 % xi=-179:2:180; yi=-89:2:90;
0309 xi = og.lon;
0310 ind = find(xi > 179.5);
0311 xi(ind) = xi(ind) - 360;
0312 yi = og.lat;
0313 [ul,vl] = uvacube2latlon(aa_regrid.x,aa_regrid.y, ...
0314 by6u,by6v, xi,yi);
0315 ave(ii ).r(iz,:,:) = ul';
0316 ave(ii+1).r(iz,:,:) = vl';
0317
0318 if flags.debug_lev > 0
0319 figure(1), subplot(1,1,1)
0320 subplot(2,1,1), surf(ul'), axis equal, view(2), shading flat
0321 subplot(2,1,2), surf(vl'), axis equal, view(2), shading flat
0322 end
0323
0324 elseif ( vars{ave(ii).ivar}{6} == 'v' )
0325
0326 % if iz == 1
0327 % disp(sprintf( ' %s was computed along with %s', ...
0328 % vars{ave(ii).ivar}{2}, vars{ave(ii-1).ivar}{2} ));
0329 % end
0330
0331 elseif ( vars{ave(ii).ivar}{6} == 'u' )
0332
0333 disp('Error: please specify regrid type for vectors');
0334 disp(' using "flags.regrid_type_for_vec"');
0335
0336 end
0337 else
0338 % conservative for scalars
0339 for jj = 1:length(rgrid.con.src_address)
0340 lld(rgrid.con.dst_address(jj)) = ...
0341 lld(rgrid.con.dst_address(jj)) ...
0342 + ( dcp(rgrid.con.src_address(jj)) ...
0343 * rgrid.con.remap_matrix(jj,1) );
0344 end
0345 ave(ii).r(iz,:,:) = reshape(lld',og.nlat,og.nlon);
0346 end
0347
0348 end
0349
0350 end
0351
0352
0353 % write the regridded fields to netCDF
0354 tcurr = 0.25 * it;
0355 nc_outf{'time'}(it-it0) = tcurr;
0356 nc_outf{'time_bnds'}(it-it0,:) = [ tcurr tcurr+0.25 ];
0357 for ii = 1:length(ave)
0358
0359 iv = ave(ii).ivar;
0360 nc_outf{ vars{iv}{2} }(it-it0,:,:) = ave(ii).r(1,:,:);
0361
0362 end
0363
0364
0365 end
0366 clear tmp
0367
0368 % close the input files
0369 for ifacet = 1:nfacets
0370 nc = close( ncl{ifacet} );
0371 end
0372 clear ncl;
0373
0374 % close output file
0375 nc_outf = close(nc_outf);
0376
0377
