% Plot average biomass distribution through the area
%

clear all
% close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inFold = '/home/martin/opendrift/PelagiaIBM/PelagiaMFS09/';
subFold = 'C2k8_1-1/';
fPrefix = 'Pelagia_Med_fwd_';
ryears = 2000:2020;
initialMass = 10000;
EUlatlim = [30 46];
EUlonlim = [-7 36.5];
saveIMG = 1;
IMGfold = 'IMG/biomass_density_plot_yearly/';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mkdir([IMGfold,subFold])

tLONS = [];
tLATS = [];
tMASS = [];

% inFiles = dir([inFold,subFold,fPrefix,'*.nc']);
n=1;
for ryear = ryears
   inFiles(n).name = [fPrefix,num2str(ryear),'.nc'];
   inFiles(n).folder = [inFold,subFold];
   n=n+1;
end


for nf = 1:length(inFiles)
    inFile = [inFiles(nf).folder,'/',inFiles(nf).name];
    disp(['Reading ',inFiles(nf).name,' ...'])
    nc = netcdf(inFile);
    lons = nc{'lon'}(:);
    lats = nc{'lat'}(:);
    nctime = nc{'time'}(1); % "seconds since 1970-01-01 00:00:00"
    biomass = nc{'biomass'}(:);
    close(nc)
    
    ryear = year(datetime(1970,1,1) + seconds(nctime));
    
    LONS = lons(:,1:1:end);
    LATS = lats(:,1:1:end);
    MASS = biomass(:,1:1:end);
    
    
    
    tLONS = max(min(LONS(:),EUlonlim(2)),EUlonlim(1));
    tLATS = max(min(LATS(:),EUlatlim(2)),EUlatlim(1));
    
    GridLat = EUlatlim(1):0.2:EUlatlim(2)+0.5;
    GridLon = EUlonlim(1)-0.5:0.2:EUlonlim(2);
    SumBiomass = zeros(length(GridLon), length(GridLat));
    
    disp('Calculating biomass distribution ...')
    for nlat = 1:length(GridLat)
        Glat = GridLat(nlat);
        for nlon = 1:length(GridLon)
            Glon = GridLon(nlon);
            idx = tLONS >= Glon-0.1 & tLONS < Glon+0.1 & tLATS >= Glat-0.1 & tLATS < Glat+0.1;
            SumBiomass(nlon,nlat) = sum(MASS(idx));
        end
    end
    
    Noutputs = size(biomass,2); % divide by the number of outputs (usually days)
    AvgBiomass = SumBiomass / Noutputs / initialMass / (0.2*0.2); % divide by number of outputs and initial mass to get average biomass per square degree

    
    disp('Plotting ...')
    hFig = figure(31);
    set(hFig, 'Position', [10 10 1200 700]); % fixed figure dimensions
    set(hFig,'Color','w')
    clf
    m_proj('mercator' ,'lon',[EUlonlim(1) EUlonlim(2)], ...
        'lat',[EUlatlim(1) EUlatlim(2)]);
    m_grid('xtick',[EUlonlim(1):4:EUlonlim(2)],'ytick',[EUlatlim(1):2:EUlatlim(2)],'linestyle','none')
    
    hold on
    m_pcolor(GridLon,GridLat,AvgBiomass')
    maxc = min(2000, max(AvgBiomass(:)));
    caxis([0 maxc])
    shading flat
    colorbar
    hold off
    m_gshhs_l()
    title(['Biomass distribution ',num2str(ryear)],'interpreter','none');
    
    
    if saveIMG
        % save figure
        fOut = [[IMGfold,subFold],'Pelagia_density_',strip(subFold,'/'),'_',num2str(nf,'%02d')];
        export_fig(fOut, '-png');
        print('-depsc','IMG/Pelagia_density_p.eps')
    end
    
    
end
