Back to home page

MITgcm

 
 

    


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