parallel.m

profile on
%% imports
file = 'o3_surface_20180701000000.nc';
ncinfo(file);
[lat] = double(ncread(file, 'lat'));
[lon] = double(ncread(file, 'lon'));
hour = ncread(file, 'hour');
bigdata= ncread(file, 'chimere_ozone');
bigdata(:,:,:,2)= ncread(file, 'emep_ozone');
bigdata(:,:,:,3)= ncread(file, 'eurad_ozone');
bigdata(:,:,:,4)= ncread(file, 'lotoseuros_ozone');
bigdata(:,:,:,5)= ncread(file, 'match_ozone');
bigdata(:,:,:,6)= ncread(file, 'mocage_ozone');
bigdata(:,:,:,7)= ncread(file, 'silam_ozone');

%% start parallel pool
NumProcessors = 4; % change this to vary the number of processors used
if isempty(gcp('nocreate')) % check if we already have a parallel pool?
    parpool(NumProcessors);
end

%% calculate mean
ozone_mean = zeros(size(bigdata,1),size(bigdata,2),size(bigdata,3));
LoopSize = size(ozone_mean,2);
LoopSize2 = size(ozone_mean,3);
parfor idx1 = 1: size(bigdata,1)
    for idx2 = 1:LoopSize
        for idx3 = 1: LoopSize2
        ozone_mean(idx1,idx2,idx3) = sum(bigdata(idx1,idx2,idx3,:))/7;
        end
    end
end

%% rotate and slice mean
o3 = zeros(size(ozone_mean,2),size(ozone_mean,1),size(ozone_mean,3));
simple_ensemble = zeros(size(ozone_mean,2)-2,size(ozone_mean,1)-2,size(ozone_mean,3));
for idx = 1:size(ozone_mean,3)
    o3(:,:,idx) = ozone_mean(:,:,idx)';
    simple_ensemble(:,:,idx)=o3(3:end,3:end,idx);
end

%% transform CBE data into o3 levels 
CBEo3 = zeros(size(CBE,1),size(CBE,2),size(CBE,3));
LoopSize = size(CBEo3,2);
LoopSize2 = size(CBEo3,3);
parfor idx1 = 1: size(CBE,1)
    for idx2 = 1:LoopSize
        for idx3 = 1: LoopSize2
        CBEo3(idx1,idx2,idx3) = CBE(idx1,idx2,idx3)*1.0497548e-07+2.9458301e-08
        end
    end
end

%% minus simple enemble from CBEo3
comparison = zeros(size(simple_ensemble,1),size(simple_ensemble,2),size(simple_ensemble,3));
LoopSize = size(comparison,2);
LoopSize2 = size(comparison,3);
parfor idx1 = 1: size(simple_ensemble,1)
    for idx2 = 1:LoopSize
        for idx3 = 1: LoopSize2
        comparison(idx1,idx2,idx3) = CBEo3(idx1,idx2,idx3)-simple_ensemble(idx1,idx2,idx3)
        end
    end
end

%% cut and transform lat and lon into the correct dimensions
lon=lon(3:end,:);
lat=lat(3:end,:);
[lon,lat] = meshgrid(lon, lat);

%% Plot contour map
for idx = 1:size(comparison,3)
    figure(idx);
    clf
    
    % create the map
    worldmap('Europe'); % set the part of the earth to show
    load coastlines
    plotm(coastlat,coastlon,'r','LineWidth',1)

    land = shaperead('landareas', 'UseGeoCoords', true);
    geoshow(gca, land, 'FaceColor', 'none','EdgeColor', 'none')

    lakes = shaperead('worldlakes', 'UseGeoCoords', true);
    geoshow(lakes, 'FaceColor', 'blue')

    rivers = shaperead('worldrivers', 'UseGeoCoords', true);
    geoshow(rivers, 'Color', 'blue')

    cities = shaperead('worldcities', 'UseGeoCoords', true);
    geoshow(cities, 'Marker', '.', 'Color', 'red')

    % display the data
    NumContours = 10;
    contourfm(lat(1:2:end,1:2:end), lon(1:2:end,1:2:end), comparison(1:2:end,1:2:end,idx) , NumContours,'LineColor','none');
    caxis([-7e-8,7e-8])
    h=colorbar;
    ylabel(h,'Difference of o3 between simple ensemble and CBE data');

    % This is a bit advanced, sets the visibility of the various parts of the
    % plot so the land, cities etc shows through.
    Plots = findobj(gca,'Type','Axes');
    Plots.SortMethod = 'depth';
    fname = sprintf ( 'Figure%i.fig', idx );
    savefig(fname);
    close all;
    fprintf('%i/25 figures complete\n', idx)
end
profile viewer