Back to home page

MITgcm

 
 

    


Warning, /verification/fizhi-cs-aqualev20/scripts/assemble_GT.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 GT fields:
                0006 %
                0007 %    For all time steps:
                0008 %      For all variables:
                0009 %    Convert units as necessary:
                0010 %    Write netCDF output with all attributes:
                0011 %
                0012 
                0013 %======================================================================
                0014 %
                0015 %  Define the connections between the APE "GT" variables and the
                0016 %  MITgcm diagnostics output
                0017 %
                0018 months.i_start = 7;
                0019 months.i_end   = 42;
                0020 bname = 'data/GTall.nc';
                0021 oname = { 'GT_fields.nc' };
                0022 
                0023 %  Get the following from:
                0024 %    ! (cd ../../ape_data_specs/ ; ./gt_parse.sh)
                0025 vars = {};
                0026 vars{1} = {'GT01','gt_sw_toai','RADSWT_ave', ...
                0027            'toa_incoming_shortwave_flux','W m-2','-','1',''};
                0028 vars{2} = {'GT02','gt_sw_toar','OSR_ave', ...
                0029            'toa_outgoing_shortwave_flux','W m-2','-','-1',''};
                0030 vars{3} = {'GT04','gt_lw_toa','OLR_ave', ...
                0031            'toa_outgoing_longwave_flux','W m-2','-','1',''};
                0032 vars{4} = {'GT07','gt_cld_frac','CLDFRC_ave', ...
                0033            'cloud_area_fraction','1','-','1',''};
                0034 vars{5} = {'GT08','gt_cldw','---', ...
                0035            'atmosphere_cloud_condensed_water_co','kg m-2','-','1',''};
                0036 vars{6} = {'GT09','gt_cldi','---', ...
                0037            'atmosphere_cloud_ice_content','kg m-2','-','1',''};
                0038 vars{7} = {'GT11','gt_cppn','PRECON_ave', ...
                0039            'convective_precipitation_flux','kg m-2 s-1','-',...
                0040            '0.000011574074','units = PRECON / (24*3600)'};
                0041 vars{8} = {'GT12','gt_dppn','GT12', ...
                0042            'large_scale_precipitation_flux','kg m-2 s-1','-', ...
                0043            '0.000011574074','GT12 = PREACC - PRECON'};
                0044 vars{9} = {'GT13','gt_evap','EVAP_ave', ...
                0045            'evaporation_flux','kg m-2 s-1','-', ...
                0046            '0.000011574074','units = EVAP / (24*3600)'};
                0047 vars{10} = {'GT15','gt_sswi','GT15', ...
                0048             'surface_downwelling_shortwave_flux','W m-2','-', ...
                0049             '1','GT15 = RADSWG/(1 - (ALBNIRDF + ALBVISDF)/2)'};
                0050 vars{11} = {'GT16','gt_sswr','GT16', ...
                0051             'surface_upwelling_shortwave_flux','W m-2','-','-1', ...
                0052             'GT16 = RADSWG*(1/(1 - (ALBNIRDF + ALBVISDF)/2) - 1)'};
                0053 vars{12} = {'GT18','gt_slwd','LWGDOWN_ave', ...
                0054             'surface_downwelling_longwave_flux','W m-2','-','-1',''};
                0055 vars{13} = {'GT19','gt_slwu','LWGUP_ave', ...
                0056             'surface_upwelling_longwave_flux','W m-2','-','1',''};
                0057 vars{14} = {'GT21','gt_slh','EFLUX_ave', ...
                0058             'surface_upward_latent_heat_flux','W m-2','-','1',''};
                0059 vars{15} = {'GT22','gt_ssh','HFLUX_ave', ...
                0060             'surface_upward_sensible_heat_flux','W m-2','-','1',''};
                0061 vars{16} = {'GT24','gt_ts2m','T2M_ave', ...
                0062             'surface_air_temperature','K','-','1',''};
                0063 vars{17} = {'GT25','gt_ps','PS_ave', ...
                0064             'surface_air_pressure','Pa','-','100',''};
                0065 
                0066 
                0067 %======================================================================
                0068 %
                0069 %  Open the input files and get the (minimum) number of time steps
                0070 %  available across all the files
                0071 %
                0072 disp('Finding available time steps')
                0073 
                0074 ncl = {};
                0075 nTmin = 0;
                0076 ifacet = 1;
                0077 % for ifacet = 1:nfacets
                0078   fname = bname;
                0079   nc = netcdf(fname, 'nowrite');
                0080   ncl{ifacet} = nc;
                0081   T = nc{'T'}(:);
                0082   nT = length(T);
                0083   if (ifacet == 1)
                0084     nTmin = nT;
                0085   else
                0086     nTmin = min(nTmin, nT);
                0087   end
                0088   nc = [];
                0089 % end
                0090 disp(sprintf('  total iterations found = %d',nTmin));
                0091 
                0092 
                0093 %======================================================================
                0094 %
                0095 %  Open the input files and get the (minimum) number of time steps
                0096 %  available across all the files
                0097 %
                0098 disp('Reading the data')
                0099 
                0100 nc = ncl{1};
                0101 time = nc{ 'T' }(:);
                0102 
                0103 for iv = 1:length(vars)
                0104       
                0105   %  cull unknown variables 
                0106   if ( vars{iv}{3}(1) == '-' )
                0107     continue
                0108   end
                0109   
                0110   %  skip over computed quantities
                0111   if strcmp( vars{iv}{1} , vars{iv}{3} )
                0112     continue
                0113   end
                0114 
                0115   %  WARN about missing variables
                0116   nc = ncl{ifacet};
                0117   if prod(size(nc{vars{iv}{3}})) == 0
                0118     if it == 1
                0119       str = 'var "%s" does not exist in file "%s"';
                0120       disp(sprintf(['  Warning: ' str], vars{iv}{3}, name(nc)));
                0121     end
                0122     continue
                0123   end
                0124       
                0125   %  get the data
                0126   t0 = squeeze( nc{ vars{iv}{3} }(:) );
                0127   
                0128   dat(iv,ifacet).n = 1;
                0129   dat(iv,ifacet).a = t0;
                0130 
                0131 end
                0132 
                0133 %  fill in the derived fields
                0134 for iv = 1:length(vars)
                0135   tmp = [];
                0136   nc = ncl{ifacet};
                0137       
                0138   if strncmp(vars{iv}{1},'GT12',4)
                0139     %  GT12 = PREACC_ave - PRECON_ave
                0140     tmp = squeeze( nc{ 'PREACC_ave' }(:) ) ...
                0141           - squeeze( nc{ 'PRECON_ave' }(:) );
                0142   end
                0143   
                0144   if  strncmp(vars{iv}{1},'GT15',4) || strncmp(vars{iv}{1},'GT16',4)
                0145     albedo = ( squeeze( nc{ 'ALBNIRDF_ave' }(:)) ...
                0146                + squeeze( nc{ 'ALBVISDF_ave' }(:)) ) ./ 2.0;
                0147     if vars{iv}{1} == 'GT15'
                0148       %  GT15 = RADSWG_ave/(1 - ALBEDO_ave)
                0149       tmp = squeeze( nc{ 'RADSWG_ave' }(:) ) ...
                0150             ./ (1.0 - albedo);
                0151     end
                0152     if strncmp(vars{iv}{1},'GT16',4)
                0153       %  GT16 = RADSWG_ave*(1/(1 - ALBEDO_ave) - 1)
                0154       tmp = squeeze( nc{ 'RADSWG_ave' }(:) ) ...
                0155             .* (1.0./(1.0 - albedo) - 1.0);
                0156     end
                0157   end
                0158   
                0159   if length(tmp) > 0
                0160     dat(iv,ifacet).n = 1;
                0161     dat(iv,ifacet).a = tmp;
                0162   end
                0163   
                0164 end
                0165 
                0166 ave = [];
                0167 for iv = 1:length(vars)
                0168   ind = length(ave) + 1;
                0169   ifacet = 1;
                0170   if length(dat(iv,ifacet).n) == 0
                0171     continue
                0172   end
                0173   ave(ind).v(:,ifacet) = dat(iv,ifacet).a;
                0174   ave(ind).ivar = iv;
                0175 end
                0176 
                0177 %  close the input files
                0178 nc = close( ncl{1} );
                0179 clear ncl;
                0180 
                0181 
                0182 %======================================================================
                0183 %
                0184 %  Convert units
                0185 %
                0186 disp('Converting units')
                0187 
                0188 for ia = 1:length(ave)
                0189   if strcmp(vars{ ave(ia).ivar }{7}, '-') ~= 1
                0190     eval(['fac = ' vars{ ave(ia).ivar }{7} ';']);
                0191     if fac ~= 1
                0192       ave(ia).v = fac .* ave(ia).v;
                0193     end
                0194   end
                0195 end
                0196 
                0197 
                0198 %======================================================================
                0199 %
                0200 %  Write netCDF output
                0201 %
                0202 disp('Writing netCDF output')
                0203 
                0204 nc = netcdf(oname{1}, 'clobber');
                0205 %nc.title = [ 'Aqua Planet: ' ...
                0206 %             'Single-Level 2-D Means from "Example" Experiment' ]; 
                0207 nc.institution = 'MIT Dept. of EAPS, Cambridge, MA, USA';
                0208 nc.source = [ 'MITgcm: ' ];
                0209 nc.Conventions = 'CF-1.0';
                0210 %nc.history = [ 'Original data produced: ' '2002/08/20' ];
                0211 
                0212 nc('time') = length( ave(1).v );
                0213 nc('bnd') = 2;
                0214 
                0215 nc{'time'} = ncdouble('time');
                0216 nc{'time'}.standard_name = 'time';
                0217 nc{'time'}.units = 'days since 0000-01-01';
                0218 nc{'time'}.bounds = 'time_bnds';
                0219 nc{'time'}(:) = time / (24*3600);
                0220 
                0221 %float time_bnds(time, bnd) ;
                0222 %time_bnds:long_name = "time interval endpoints" ;
                0223 nc{'time_bnds'} = ncdouble('time','bnd');
                0224 nc{'time_bnds'}.ape_name = 'time interval endpoints';
                0225 nc{'time_bnds'}.units = 'days since 0000-01-01';
                0226 tbds(1,:) = (time - 24*3600)/(24*3600);
                0227 tbds(2,:) = (time)/(24*3600);
                0228 nc{'time_bnds'}(:) = tbds';
                0229 
                0230 
                0231 for ii = 1:length(ave)
                0232   
                0233   iv = ave(ii).ivar;
                0234   nc{ vars{iv}{2} } = ncfloat( 'time' );
                0235   nc{ vars{iv}{2} }.ape_name = vars{iv}{4};
                0236   nc{ vars{iv}{2} }.units = vars{iv}{5};
                0237   nc{ vars{iv}{2} }.FillValue_ = 1.0e20;
                0238   
                0239   nc{ vars{iv}{2} }(:) = ave(ii).v;
                0240 
                0241 end
                0242 
                0243 nc = close(nc);