Warning, /verification/fizhi-cs-aqualev20/scripts/assemble_SH.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 SH01 -- SH30 fields:
0006 %
0007 % For all time steps:
0008 % For all variables:
0009 % For all facets:
0010 % compute the "derived" fields:
0011 % sum:
0012 % Compute the averages from the sums:
0013 % Rotate vector components as necessary:
0014 % Convert units as necessary:
0015 % Regrid onto lat-lon:
0016 % Write netCDF output with all attributes:
0017 %
0018
0019 %======================================================================
0020 %
0021 % Define the connections between the APE "SH" variables and the
0022 % MITgcm diagnostics output
0023 %
0024 nfacets = 6;
0025 months.i_start = 7;
0026 months.i_end = 42;
0027 bname = 'data/SHave.t%03d.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 = { 'SH_fields.nc' };
0032
0033 flags.regrid_type_for_vec = 2;
0034 flags.debug = 0;
0035
0036 % Get the following from:
0037 % ! (cd ../../ape_data_specs/ ; ./sh_parse.sh)
0038 vars = {};
0039 vars{1} = {'SH01','sh_sw_toai','RADSWT','toa_incoming_shortwave_flux','W m-2','-','1'};
0040 vars{2} = {'SH02','sh_sw_toar','OSR','toa_outgoing_shortwave_flux','W m-2','-','-1'};
0041 vars{3} = {'SH04','sh_lw_toa','OLR','toa_outgoing_longwave_flux','W m-2','-','1'};
0042 vars{4} = {'SH07','sh_cld_frac','CLDFRC','cloud_area_fraction','1','-','1'};
0043 vars{5} = {'SH08','sh_cldw','---','atmosphere_cloud_condensed_water_co','kg m-2','-','1'};
0044 vars{6} = {'SH09','sh_cldi','---','atmosphere_cloud_ice_content','kg m-2','-','1'};
0045 vars{7} = {'SH11','sh_cppn','PRECON','convective_precipitation_flux','kg m-2 s-1','-','0.000011574074'};
0046 vars{8} = {'SH12','sh_dppn','SH12','large_scale_precipitation_flux','kg m-2 s-1','-','0.000011574074'};
0047 vars{9} = {'SH13','sh_evap','EVAP','evaporation_flux','kg m-2 s-1','-','0.000011574074'};
0048 vars{10} = {'SH15','sh_sswi','SH15','surface_downwelling_shortwave_flux','W m-2','-','1'};
0049 vars{11} = {'SH16','sh_sswr','SH16','surface_upwelling_shortwave_flux','W m-2','-','-1'};
0050 vars{12} = {'SH18','sh_slwd','LWGDOWN','surface_downwelling_longwave_flux','W m-2','-','-1'};
0051 vars{13} = {'SH19','sh_slwu','LWGUP','surface_upwelling_longwave_flux','W m-2','-','1'};
0052 vars{14} = {'SH21','sh_slh','EFLUX','surface_upward_latent_heat_flux','W m-2','-','1'};
0053 vars{15} = {'SH22','sh_ssh','HFLUX','surface_upward_sensible_heat_flux','W m-2','-','1'};
0054 vars{16} = {'SH24','sh_ts2m','T2M','surface_air_temperature','K','-','1'};
0055 vars{17} = {'SH25','sh_q2m','Q2M','specific_humidity','1','-','0.001'};
0056 vars{18} = {'SH26','sh_tauu','UFLUX','surface_downward_eastward_stress','Pa','u','-1'};
0057 vars{19} = {'SH27','sh_tauv','VFLUX','surface_downward_northward_stress','Pa','v','-1'};
0058 vars{20} = {'SH28','sh_u10m','U10M','eastward_wind','m s-1','u','1'};
0059 vars{21} = {'SH29','sh_v10m','V10M','northward_wind','m s-1','v','1'};
0060 vars{22} = {'SH30','sh_ps','PS','surface_air_pressure','Pa','-','100'};
0061
0062
0063
0064 %======================================================================
0065 %
0066 % Open the input files and get the (minimum) number of time steps
0067 % available across all the files
0068 %
0069 disp('Finding available time steps')
0070
0071 ncl = {};
0072 nTmin = 0;
0073 ifacet = 1;
0074 for ifacet = 1:nfacets
0075 fname = sprintf(bname,ifacet);
0076 nc = netcdf(fname, 'nowrite');
0077 ncl{ifacet} = nc;
0078 T = nc{'T'}(:);
0079 nT = length(T);
0080 if (ifacet == 1)
0081 nTmin = nT;
0082 else
0083 nTmin = min(nTmin,nT);
0084 end
0085 nc = [];
0086 end
0087 disp(sprintf(' total iterations found = %d',nTmin));
0088
0089
0090 %======================================================================
0091 %
0092 % Compute the derived fields and the temporal sums of all fields on
0093 % the original MITgcm grid
0094 %
0095 disp('Computing Sums')
0096
0097 dat = [];
0098 iz = 1;
0099 it = 1;
0100 iv = 1;
0101 ifacet = 1;
0102 for it = [ 1 ] % months.i_start:min(nTmin,months.i_end)
0103 disp(sprintf(' iteration = %d',it));
0104 for iv = 1:length(vars)
0105 for ifacet = 1:nfacets
0106
0107 % cull unknown variables
0108 if ( vars{iv}{3}(1) == '-' )
0109 continue
0110 end
0111 % disp(sprintf('%d, %d',iv,ifacet));
0112
0113 % skip over computed quantities
0114 if strcmp( vars{iv}{1} , vars{iv}{3} )
0115 continue
0116 end
0117
0118 % WARN about missing variables
0119 nc = ncl{ifacet};
0120 if prod(size(nc{vars{iv}{3}})) == 0
0121 if it == months.i_start
0122 str = 'var "%s" does not exist in file "%s"';
0123 disp(sprintf([' Warning: ' str], vars{iv}{3}, name(nc)));
0124 end
0125 continue
0126 end
0127
0128 if it == 1
0129 dat(iv,ifacet).n = 1;
0130 dat(iv,ifacet).a = squeeze( nc{ vars{iv}{3} }(it,iz,:,:) );
0131 else
0132 dat(iv,ifacet).n = dat(iv,ifacet).n + 1;
0133 dat(iv,ifacet).a = dat(iv,ifacet).a + ...
0134 squeeze( nc{ vars{iv}{3} }(it,iz,:,:) );
0135 end
0136
0137 end
0138 end
0139
0140 % fill in the derived fields
0141 for iv = 1:length(vars)
0142 for ifacet = 1:nfacets
0143 tmp = [];
0144 nc = ncl{ifacet};
0145
0146 if strncmp(vars{iv}{1},'SH12',4)
0147 % SH12 = PREACC - PRECON
0148 tmp = squeeze( nc{ 'PREACC' }(it,iz,:,:) ) ...
0149 - squeeze( nc{ 'PRECON' }(it,iz,:,:) );
0150 end
0151
0152 if strncmp(vars{iv}{1},'SH15',4) || strncmp(vars{iv}{1},'SH16',4)
0153 albedo = ( squeeze( nc{ 'ALBNIRDF' }(it,iz,:,:)) ...
0154 + squeeze( nc{ 'ALBVISDF' }(it,iz,:,:)) ) ./ 2.0;
0155 if vars{iv}{1} == 'SH15'
0156 % SH15 = RADSWG/(1 - ALBEDO)
0157 tmp = squeeze( nc{ 'RADSWG' }(it,iz,:,:) ) ...
0158 ./ (1.0 - albedo);
0159 end
0160 if strncmp(vars{iv}{1},'SH16',4)
0161 % SH16 = RADSWG*(1/(1 - ALBEDO) - 1)
0162 tmp = squeeze( nc{ 'RADSWG' }(it,iz,:,:) ) ...
0163 .* (1.0./(1.0 - albedo) - 1.0);
0164 end
0165 end
0166
0167 if prod(size(tmp)) > 0
0168 if it == 1
0169 dat(iv,ifacet).n = 1;
0170 dat(iv,ifacet).a = tmp;
0171 else
0172 dat(iv,ifacet).n = dat(iv,ifacet).n + 1;
0173 dat(iv,ifacet).a = dat(iv,ifacet).a + tmp;
0174 end
0175 end
0176
0177 end
0178 end
0179 end
0180 clear tmp
0181
0182 % close the input files
0183 for ifacet = 1:nfacets
0184 nc = close( ncl{ifacet} );
0185 end
0186 clear ncl;
0187
0188 %======================================================================
0189 %
0190 % Compute the averages from the temporal sums
0191 %
0192 disp('Computing Averages from Sums')
0193
0194 ave = [];
0195 for iv = 1:length(vars)
0196 ind = length(ave) + 1;
0197 for ifacet = 1:nfacets
0198 if length(dat(iv,ifacet).n) == 0
0199 continue
0200 end
0201 ave(ind).v(:,:,ifacet) = dat(iv,ifacet).a ./ dat(iv,ifacet).n;
0202 ave(ind).ivar = iv;
0203 end
0204 end
0205
0206 clear dat;
0207
0208 %======================================================================
0209 %
0210 % Rotate vector components & convert units
0211 %
0212 clear ig;
0213 gvars = { 'XC','YC','dxF','dyF','rA','XG','YG','dxV', ...
0214 'dyU','rAz','dxC','dyC','rAw','rAs','dxG','dyG' };
0215 for iface = 1:6
0216 fname = sprintf(gname, iface);
0217 nc = netcdf(fname, 'nowrite');
0218 if iface == 1
0219 ig.ne = length(nc('X'));
0220 end
0221 for jj = 1:length(gvars)
0222 comm = [ 'ig.' gvars{jj} '(:,:,iface) = nc{''' ...
0223 gvars{jj} '''}(:);'];
0224 eval(comm);
0225 end
0226 nc = close(nc);
0227 end
0228
0229 ne = ig.ne;
0230 n1 = ne - 1;
0231 dux = zeros(size(ig.XC));
0232 duy = zeros(size(ig.XC));
0233 dvx = zeros(size(ig.XC));
0234 dvy = zeros(size(ig.XC));
0235 % dux(:,:,:) = diff(ig.XG(:,1:ne,:),1,1);
0236 % dvx(:,:,:) = diff(ig.XG(1:ne,:,:),1,2);
0237 % duy(:,:,:) = diff(ig.YG(:,1:ne,:),1,1);
0238 % dvy(:,:,:) = diff(ig.YG(1:ne,:,:),1,2);
0239 %
0240 % Note the first two dimensions below are permuted relative to the
0241 % ordering in the *.mitgrid files (commented-out above).
0242 dux(:,:,:) = diff(ig.XG(1:ne,:,:),1,2);
0243 dvx(:,:,:) = diff(ig.XG(:,1:ne,:),1,1);
0244 duy(:,:,:) = diff(ig.YG(1:ne,:,:),1,2);
0245 dvy(:,:,:) = diff(ig.YG(:,1:ne,:),1,1);
0246 dux = dux + 360*double(dux < 180);
0247 dux = dux - 360*double(dux > 180); % squeeze([ min(min(dux)) max(max(dux)) ])
0248 duy = duy + 360*double(duy < 180);
0249 duy = duy - 360*double(duy > 180); % squeeze([ min(min(duy)) max(max(duy)) ])
0250 dvx = dvx + 360*double(dvx < 180);
0251 dvx = dvx - 360*double(dvx > 180); % squeeze([ min(min(dvx)) max(max(dvx)) ])
0252 dvy = dvy + 360*double(dvy < 180);
0253 dvy = dvy - 360*double(dvy > 180); % squeeze([ min(min(dvy)) max(max(dvy)) ])
0254 dut = sqrt(dux.^2 + duy.^2); % squeeze([ min(min(dut)) max(max(dut)) ])
0255 dvt = sqrt(dvx.^2 + dvy.^2); % squeeze([ min(min(dvt)) max(max(dvt)) ])
0256 ig.llux = dux ./ dut; % squeeze([ min(min(ig.llux)) max(max(ig.llux)) ])
0257 ig.lluy = duy ./ dut; % squeeze([ min(min(ig.lluy)) max(max(ig.lluy)) ])
0258 ig.llvx = dvx ./ dvt; % squeeze([ min(min(ig.llvx)) max(max(ig.llvx)) ])
0259 ig.llvy = dvy ./ dvt; % squeeze([ min(min(ig.llvy)) max(max(ig.llvy)) ])
0260 clear dux duy dvx dvy dut dvt ne n1
0261
0262 if flags.regrid_type_for_vec ~= 2
0263 % Vector rotation
0264 for ia = 1:length(ave)
0265 if strcmp(vars{ ave(ia).ivar }{6}, 'u')
0266 % llu = u .* llux + v .* llvx;
0267 % llv = u .* lluy + v .* llvy;
0268 llu = ave(ia).v .* ig.llux + ave(ia+1).v .* ig.llvx;
0269 llv = ave(ia).v .* ig.lluy + ave(ia+1).v .* ig.llvy;
0270 ave(ia).v = llu;
0271 ave(ia+1).v = llv;
0272 end
0273 end
0274 end
0275 clear llu llv
0276
0277 % Units conversion
0278 for ia = 1:length(ave)
0279 if strcmp(vars{ ave(ia).ivar }{7}, '-') ~= 1
0280 eval(['fac = ' vars{ ave(ia).ivar }{7} ';']);
0281 if fac ~= 1
0282 ave(ia).v = fac .* ave(ia).v;
0283 end
0284 end
0285 end
0286
0287 %======================================================================
0288 %
0289 % Regrid
0290 %
0291 % JMC suggested an evenly spaced Lat-Lon at: 128x64
0292 %
0293 disp('Regridding')
0294
0295 og = [];
0296 og.nlat = 64;
0297 og.nlon = 128;
0298 og.latcell = 180/og.nlat;
0299 og.loncell = 360/og.nlon;
0300 og.lat = linspace(-90+og.latcell/2, 90-og.latcell/2, og.nlat);
0301 og.lon = linspace(0+og.loncell/2, 360-og.loncell/2, og.nlon);
0302 og.latbnd(:,1) = og.lat - og.latcell/2.0;
0303 og.latbnd(:,2) = og.lat + og.latcell/2.0;
0304 og.lonbnd(:,1) = og.lon - og.loncell/2.0;
0305 og.lonbnd(:,2) = og.lon + og.loncell/2.0;
0306
0307 rgvar = { 'src_address', 'dst_address', 'remap_matrix' };
0308 rtype = { { 'con', r_con_name}, {'dwt', r_dwt_name } };
0309 for kk = 1:length(rtype)
0310 nc = netcdf(rtype{kk}{2}, 'nowrite');
0311 for jj = 1:length(rgvar)
0312 comm = sprintf('rgrid.%s.%s = nc{''%s''}(:);',...
0313 rtype{kk}{1},rgvar{jj},rgvar{jj});
0314 eval(comm);
0315 % disp(comm);
0316 end
0317 nc = close(nc);
0318 end
0319
0320 if ( flags.regrid_type_for_vec > 0 )
0321 % >> x=rdmds('XC');
0322 % >> y=rdmds('YC');
0323 % >> t=rdmds('Ttave.0000513360');
0324 % >> xi=-179:2:180;yi=-89:2:90;
0325 % >> del=cube2latlon_preprocess(x,y,xi,yi);
0326 % >> ti=cube2latlon_fast(del,t);
0327 aa_regrid.x = rdmds('aa_regrid/XC');
0328 aa_regrid.y = rdmds('aa_regrid/YC');
0329 % u = rdmds('aa_regrid/U.0000000000');
0330 % v = rdmds('aa_regrid/V.0000000000');
0331 % del = cube2latlon_preprocess(LON,LAT,xc,yc);
0332 %
0333 % aa_regrid.del = ...
0334 % cube2latlon_preprocess(aa_regrid.x,aa_regrid.y, ...
0335 % og.lon,og.lat);
0336 end
0337
0338 % save zzz
0339 % zzz
0340
0341 ii = 1;
0342 for ii = 1:length(ave)
0343
0344 dcp = zeros(prod(size( ave(ii).v )),1);
0345 nn = 0;
0346 for kk = 1:size( ave(ii).v, 3 )
0347 for jk = 1:size( ave(ii).v, 2 )
0348 for ik = 1:size( ave(ii).v, 1 )
0349 nn = nn + 1;
0350 dcp(nn) = ave(ii).v(ik,jk,kk);
0351 end
0352 end
0353 end
0354 lld = zeros(og.nlat*og.nlon,1);
0355
0356 if ( ( flags.regrid_type_for_vec > 0 ) ...
0357 && ( vars{ave(ii).ivar}{6} == 'u' ...
0358 || vars{ave(ii).ivar}{6} == 'v' ) )
0359
0360 if flags.regrid_type_for_vec == 1
0361
0362 % use distributed weights for vectors
0363 for jj = 1:length(rgrid.dwt.src_address)
0364 lld(rgrid.dwt.dst_address(jj)) = ...
0365 lld(rgrid.dwt.dst_address(jj)) ...
0366 + ( dcp(rgrid.dwt.src_address(jj)) ...
0367 * rgrid.dwt.remap_matrix(jj,1) );
0368 end
0369 ave(ii).r = reshape(lld',og.nlat,og.nlon);
0370
0371 elseif ( flags.regrid_type_for_vec == 2 ...
0372 && vars{ave(ii).ivar}{6} == 'u' )
0373
0374 % use Alistair's uvcube2latlon() for vectors
0375 nsiz = size(ave(ii).v);
0376 by6u = zeros([ nsiz(1)*nsiz(3) nsiz(2) ]);
0377 by6v = zeros([ nsiz(1)*nsiz(3) nsiz(2) ]);
0378 for kk = 1:size( ave(ii).v, 3 )
0379 is = 1 + (kk-1)*nsiz(1);
0380 ie = is + nsiz(1) - 1;
0381 by6u(is:ie,:) = ave(ii ).v(:,:,kk)';
0382 by6v(is:ie,:) = ave(ii+1).v(:,:,kk)';
0383 end
0384 % xi=-179:2:180; yi=-89:2:90;
0385 xi = og.lon;
0386 ind = find(xi > 179.5);
0387 xi(ind) = xi(ind) - 360;
0388 yi = og.lat;
0389 [ul,vl] = uvacube2latlon(aa_regrid.x,aa_regrid.y, ...
0390 by6u,by6v, xi,yi);
0391 ave(ii ).r = ul';
0392 ave(ii+1).r = vl';
0393
0394 if flags.debug > 0
0395 figure(1), subplot(1,1,1)
0396 subplot(2,1,1), surf(ul'), axis equal, view(2), shading flat
0397 subplot(2,1,2), surf(vl'), axis equal, view(2), shading flat
0398 end
0399
0400 elseif ( vars{ave(ii).ivar}{6} == 'v' )
0401
0402 disp(sprintf( ' %s was computed along with %s', ...
0403 vars{ave(ii).ivar}{2}, vars{ave(ii-1).ivar}{2} ));
0404
0405 elseif ( vars{ave(ii).ivar}{6} == 'u' )
0406
0407 disp('Error: please specify regrid type for vectors');
0408 disp(' using "flags.regrid_type_for_vec"');
0409
0410 end
0411 else
0412 % conservative for scalars
0413 for jj = 1:length(rgrid.con.src_address)
0414 lld(rgrid.con.dst_address(jj)) = ...
0415 lld(rgrid.con.dst_address(jj)) ...
0416 + ( dcp(rgrid.con.src_address(jj)) ...
0417 * rgrid.con.remap_matrix(jj,1) );
0418 end
0419 ave(ii).r = reshape(lld',og.nlat,og.nlon);
0420 end
0421
0422 if flags.debug > 0
0423 figure(2)
0424 surf(og.lon,og.lat,ll128x64), view(2), shading flat
0425 end
0426
0427 end
0428
0429
0430 %======================================================================
0431 %
0432 % Write netCDF output
0433 %
0434 disp('Writing netCDF output for SH')
0435
0436 nc = netcdf(oname{1}, 'clobber');
0437 nc.title = [ 'Aqua Planet: ' ...
0438 'Single-Level 2-D Means from "Example" Experiment' ];
0439 nc.institution = 'MIT Dept. of EAPS, Cambridge, MA, USA';
0440 nc.source = [ 'MITgcm: ' ];
0441 nc.Conventions = 'CF-1.0';
0442 nc.history = [ 'Original data produced: ' '2002/08/20' ];
0443
0444 nc('lon') = og.nlon;
0445 nc('lat') = og.nlat;
0446 nc('bnd') = 2;
0447
0448 nc{'lon'} = ncdouble('lon');
0449 nc{'lon'}.standard_name = 'longitude';
0450 nc{'lon'}.units = 'degrees_east';
0451 nc{'lon'}.bounds = 'lon_bnds';
0452 nc{'lon'}(:) = og.lon;
0453
0454 nc{'lat'} = ncdouble('lat');
0455 nc{'lat'}.standard_name = 'latitude';
0456 nc{'lat'}.units = 'degrees_north';
0457 nc{'lat'}.bounds = 'lat_bnds';
0458 nc{'lat'}(:) = og.lat;
0459
0460 nc{'lon_bnds'} = ncdouble('lon','bnd');
0461 nc{'lon_bnds'}.ape_name = 'longitude cell bounds';
0462 nc{'lon_bnds'}.units = 'degrees_east';
0463 nc{'lon_bnds'}(:) = og.lonbnd;
0464
0465 nc{'lat_bnds'} = ncdouble('lat','bnd');
0466 nc{'lat_bnds'}.ape_name = 'latitude cell bounds';
0467 nc{'lat_bnds'}.units = 'degrees_north';
0468 nc{'lat_bnds'}(:) = og.latbnd;
0469
0470 for ii = 1:length(ave)
0471
0472 iv = ave(ii).ivar;
0473 nc{ vars{iv}{2} } = ncfloat( 'lat', 'lon' );
0474 nc{ vars{iv}{2} }.ape_name = vars{iv}{4};
0475 nc{ vars{iv}{2} }.units = vars{iv}{5};
0476 nc{ vars{iv}{2} }.FillValue_ = 1.0e20;
0477
0478 % deal with missing values
0479 ind = find(isnan( ave(ii).r ));
0480 ave(ii).r(ind) = 1.0e20;
0481
0482 nc{ vars{iv}{2} }(:) = ave(ii).r;
0483
0484 end
0485
0486 nc = close(nc);
0487
