function [blocking, time, eval_lon, eval_lat] = blocking_frequency_theta(path,mon,y, filt_persist, large_scale, noleap)
% Function to calculate 2D blocking index, following Woollings et al. (2008).
% Uses loadnc.m, which simply opens a netcdf file
%       and selects time steps / periods / months
% Uses weight_matrix, that computes the relative size of each grid cell
%       and creates a weighting array, to be multiplied with the data
% G. Gollan, 2017, ggollan@geomar.de

if ~exist('filt_persist','var')
    filt_persist = 1;
end

if ~exist('large_scale','var')
    large_scale = 1;
end

% throw out leap days? 
if ~exist('noleap','var')
    noleap = 0;
end

% load data
[theta,lon,lat,time] = loadnc(path,0,mon,y, 'noleap', noleap); % ,'nol',-1); % ,'center',0);
[~,~,T] = size(theta);

% central latitudes and longitudes at which the meridional gradient shall be computed
eval_lat = (25 : 4 : 73);
clatte = ones( length(lon) ,1);
clatte = clatte * eval_lat;

eval_lon = lon(1) : 5 : lon(end);

% pre-allocate arrays
x = length( eval_lon );
y = length( eval_lat );

blocking      = zeros(x, y, T);
blocking_inst = zeros(x, y, T);
theta_tmp     = zeros(x, length(lat), T);

% longitude ranges to compute zonal mean. Works for lon-arrays centred around 0 or 180 deg
for i = 1 : x
    clon = eval_lon(i);
    if clon <= - 177.5
        sellon = find(lon >= clon + 360 - 2.5 | lon <= clon + 2.5 );
    elseif clon >= 177.5
        sellon = find(lon >= clon - 2.5 | lon <= clon - 360 + 2.5 );
    else
        sellon = find(lon >= clon - 2.5 & lon <= clon + 2.5 );
    end
    theta_tmp(i, :, :) = squeeze(mean( theta(sellon, :, :), 1 ));
end

% weight data according to cell size (in case of lon-lat grids, this is the cosine of the latitude)
% to compute area-averages
w_tmp = weight_matrix( zeros(x, length(lat)), lat, 0);
theta_tmp = bsxfun(@times, w_tmp, theta_tmp);

w = weight_matrix( theta(:,:,1), lat, 0);
for j = 1: y
    if j == 3
        'gh';
    end
    clat = clatte(:,j);
    d_phi = 15;
    phi_n = clat + d_phi;
    phi_0 = clat        ;
    phi_s = clat - d_phi;

    south = zeros( x, T );
    north = south;
    for i = 1 : x

        south(i, :) = squeeze( sum(     theta_tmp( i, lat >= phi_s(i) & lat <= phi_0(i), :), 2 ) ); % 1 ./ d_phi .* 
        south(i, :) = south(i, :) ./ nansum( w( 1, lat >= phi_s(i) & lat <= phi_0(i)   ) );

        north(i, :) = squeeze( sum(     theta_tmp( i, lat >= phi_0(i) & lat <= phi_n(i), :), 2 ) ); % 1 ./ d_phi .*
        north(i, :) = north(i, :) ./ nansum( w( 1, lat >= phi_0(i) & lat <= phi_n(i)   ) );

    end

    blocking_tmp = zeros(x, T);
    blocking_tmp(north - south > 0) = 1;
    
    blocking(:, j, :) = blocking_tmp;
    
    inst = north - south;
    blocking_inst(:, j, :) = inst;
end

% apply spatial and temporal contraints as in Woollings et al. 2008
if large_scale
    % take maximum of instantaneous index within +-4° of latitude
    blck_tmp = zeros(size( blocking_inst ));

    for j = 1:y
        if j == 1
            blck_tmp(:,j,:) = max(blocking_inst(:, j: j+1, :), [], 2);
        elseif j == y
            blck_tmp(:,j,:) = max(blocking_inst(:, j-1 :j, :), [], 2);
        else
            blck_tmp(:,j,:) = max(blocking_inst(:, j-1:j+1, :), [], 2);
        end
    end

    
    % extract "large scale wave breaking" (at least 15° of longitude [3 adjacent longitudes])
    blck_tmp = cat(1, blck_tmp(end,:,:), blck_tmp, blck_tmp(1,:,:) ); % merge longitudes
    blck_large = zeros(size(blck_tmp));
    x = length(blck_tmp(:,1,1));
    
    for i = 2:x-1
        for j = 1:y
            for t = 1:T
                if  blck_tmp  (i-1 : i+1, j, t) > 0
                    blck_large(i-1 : i+1, j, t) = 1;
                end
            end
        end
    end
    blocking = blck_large( 2:end-1, :, :);
end

% W08: a wave-breaking episode is identified at a given point if
% large-scale wave breaking occurs within 10° of longitude of that point
% for five consecutive days.
episode = 1;
if episode && large_scale
    filt_persist = 1;
    n = 5;
    blck_tmp = cat(1, blocking(end-1:end,:,:), blocking, blocking(1:2,:,:) ); % merge longitudes
    blck_episode = zeros(size(blck_tmp)); 
    x = length(blck_tmp(:,1,1));
    
    for i = 3:x-2
        for j = 1:y
            for t = 1:T
                % to avoid trying to access array points outside of
                % range
                if t+n > T
                    break
                end  
                tmp = squeeze( blck_tmp(i-2:i+2, j, t:t+n-1) );
                if  sum( tmp, 1 ) >= 1
                    blck_episode(i, j, t) = 1;
                end
            end
        end
    end
    blocking_epi = blck_episode( 3:end-2, :, :);
end

%%%%%%%%%%%%%%%
% blocking episode 
%%%%%%%%%%%%%%%
blocking = blocking_epi;

end