Back to home page

MITgcm

 
 

    


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