function [ emissionsEDGAR, latlon ] = prepareEDGAR( edgarDir, topcorner, botcorner, ...
                                            plotFigs )
%prepareEDGAR Reads, crops and projects EDGAR emissions to specified region
%   emissionsEDGAR = prepareEDGAR( edgarDir, tc, bc ) reads in areas of
%   cells, their locations and the actual emissions in 1deg x 1deg cells
%   (EDGAR), worldwide, and then crop out the region specified by the
%   topcorner and botcorner lat/lon in tc, bc
%
% INPUTS :
%   edgarDir : data directory with area_1x1.nc, and other inputs
%   tc       : We want edgar emissions in a rectangular region, broken into
%              1 deg x 1 deg boxes. tc contains the lat/lon of the cell
%              center of the top left box. lats run -90 to 90 and longs
%              -180 to 180
%   bc       : ditto, but bottom right box. 
%   plotFigs : true or false
%
% OUTPUTS :
%   emissionsEDGAR : a 2D array, containing CO2 emissions at 1 deg x 1deg 
%                    cell-centers in the region shown by tc, bc. The units
%                    micromoles co2/sec/m^2
%   latlon : 3D array. latlon(:, :, 1) contains the lats of the
%            cell-centers of each cell, latlon(:, :, 2) the longitudes
%
% Jaideep Ray, 7/12/2012
% -------------------------------------------------------------------------

% Edgar stores emissions as (-89.5:89.5, 0.5:359.5), and areas go
% (-89.5:89.5, -179.5:-179.5)
nlong = 360 ; % Cells in the longitudinal direction; counting east from Greenwich
nlat  = 180 ; % Cells in the latitudinal direction; counting north from South Pole
latfile  = [edgarDir, '/latitude_1x1.dat'] ;
lonfile  = [edgarDir, '/longitude_1x1.dat'] ;
areafile = [edgarDir, '/cell_area_1x1.dat'] ;
emifile  = [edgarDir, '/emi_co2_1x1.dat'] ;
aviafile = [edgarDir, '/emi_co2_aviation_1x1.dat'] ;

%% ---- Read the files
[h, d] = hdrload(latfile)  ; tmp = d' ; lat   = tmp(:) ; % latitudes of cell-centers (N)
[h, d] = hdrload(lonfile)  ; tmp = d' ; lon   = tmp(:) ; % longitudes of cell-centers (M)
[h, d] = hdrload(areafile) ; tmp = d' ; areas = tmp(:) ; % areas of cell-centers (N x M)
[h, d] = hdrload(emifile)  ; tmp = d' ; emi   = tmp(:) ; % emissions (kg/m^2/s) in each cell (N x M)
[h, d] = hdrload(aviafile) ; tmp = d' ; emia  = tmp(:) ; % airplane emissions (kg/m^2/s) in each cell (N x M)

%% ---- Turn the data into 2D arrays
areas2D = reshape(areas, nlong, nlat)  ; areas2D  = transpose(areas2D) ; % 180 rows X 360 cols
emi2D   = reshape(emi,   nlong, nlat)  ; emi2D    = transpose(emi2D)   ; % 180 x 360 cols
emia2D  = reshape(emia,   nlong, nlat) ; emia2D   = transpose(emia2D)   ; % 180 x 360 cols

% Fix the areas file. While emissions are stored (-89.5:89.5, 0.5:359.5),
% areas go (-89.5:89.5, -179.5:-179.5). So the first 180 cols of areas
% should be moved to the back
tmp = areas2D(:, 1:end/2) ;
areas2D(:, 1:end/2) = areas2D(:, end/2+1:end) ;
areas2D(:, end/2+1:end) = tmp ;

%% %% ---- Plot to see if OK. Remove aviation fluxes

emissions = areas2D .* (emi2D - emia2D) / 1000 ; % Tonnes of co2 per grid-cell per second
emissions = emissions * 3600    ; % tonnes of co2 per grid-cell per hour
ii = emissions < 1.0e-10 ; emissions(ii) = 1.0e-10 ;
log10emissions = log10(emissions) ;

[longs_grid, lats_grid] = meshgrid(lon, lat) ;
if (plotFigs == true)
    figure(1) ;
    pcolor(longs_grid, lats_grid, log10emissions) ; shading flat ; colorbar ;
    caxis([0, 0.8*max(log10emissions(:))]);
    title('log_{10}(emissions) [tonnes/hr/gridcell] (EDGAR, 2005, annually averaged)') ;
    xlabel('Longitude, 1 degree resolution') ;
    ylabel('Latitude, 1 degree resolution') ;
end

totalEmissionsEDGAR = sum(emissions(:)) * 24 * 365 ;
fprintf(' ===> Total **world** emissions in 2005, from EDGAR = %e tonnes of CO2\n', ...
        totalEmissionsEDGAR);

%% Extract EDGAR fluxes over the domain favored by vulcan

% Turn Gurney's numbering scheme into one favored by EDGAR
tc = topcorner ; tc(2) = 360 + tc(2) ;
bc = botcorner ; bc(2) = 360 + bc(2) ;
lat_gurney  = bc(1):1.0:tc(1) ;
long_gurney = tc(2):1.0:bc(2) ;
[longs_grid_gurney, lats_grid_gurney] = meshgrid(long_gurney, lat_gurney) ;
log10US_emissions = interp2(longs_grid, lats_grid, log10emissions, ...
                            longs_grid_gurney, lats_grid_gurney) ;
US_emissions = 10.^log10US_emissions ;         

if (plotFigs == true)
    figure(2) ;
    pcolor(longs_grid_gurney, lats_grid_gurney, log10US_emissions) ; shading flat ; colorbar ;
    caxis([0, 0.8*max(log10emissions(:))]);
    title('log_{10}(emissions) [tonnes/hr/gridcell] (EDGAR, 2005, annually averaged)') ;
    xlabel('Longitude, 1 degree resolution') ;
    ylabel('Latitude, 1 degree resolution') ;
end

totalEmissionsVulcan = sum(US_emissions(:)) * 24 * 365 ;
fprintf(' ===> Total **US** emissions in 2005, from EDGAR = %e tonnes of CO2\n', ...
        totalEmissionsVulcan);

%% ---- Now, redo, but in terms of fluxes
emissions = (emi2D - emia2D) ; % kgs co2 per grid-cell per second ;
ii = emissions < 1.0e-12 ; emissions(ii) = 1.0e-20 ;
log10emissions = log10(emissions) ;

% Turn Gurney's numbering scheme into one favored by EDGAR
tc = topcorner ; tc(2) = 360 + tc(2) ;
bc = botcorner ; bc(2) = 360 + bc(2) ;
lat_gurney  = bc(1):1.0:tc(1) ;
long_gurney = tc(2):1.0:bc(2) ;
[longs_grid_gurney, lats_grid_gurney] = meshgrid(long_gurney, lat_gurney) ;
log10US_emissions = interp2(longs_grid, lats_grid, log10emissions, ...
                            longs_grid_gurney, lats_grid_gurney) ;
US_emissions = 10.^log10US_emissions ;        

% ---- Now convert back to micromoles/m^2/sec
MolarWeightCO2 = (12+32) * 1.0e-3 ; % 0.044 kg/mole for CO2
microMolarWeight = MolarWeightCO2 / 1.0e6 ; % a million micromoles in 1 molar wight (kg/mu_Mole)
emissionsEDGAR = US_emissions / microMolarWeight ;

latlon(:, :, 1) = lats_grid_gurney ;
latlon(:, :, 2) = longs_grid_gurney ;
end

