function [out]=PR(dx,H,f,Nm)
% Compute modes using Proudman-Rao method. Compute rotating and non
% rotating modes as basis for generalized modes.

ii=complex(0,1);
g=9.81;
[Nx,Ny]=size(H);

% Find minimum nonzero depth
keep=H>0;
H_min=min(H(keep));

% Find H at grid corners for psi
H2=zeros([Nx-1,Ny-1]);
for i=1:Nx-1
    for j=1:Ny-1
        if H(i,j)~=0 && H(i+1,j)~=0 && H(i,j+1)~=0 && H(i+1,j+1)~=0
            H2(i,j)=(H(i,j)+H(i+1,j)+H(i,j+1)+H(i+1,j+1))/4;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define L1 matrix div(H grad()) with div(phi).N=0 on boundary
water=find(H>0); % Linear index to reconstruct modes
[coord.i,coord.j]=find(H>0); % Row/column indicies to build matrix
N=length(coord.i);

L1=sparse(N,N); 
for ind=1:N
    
    % The coordinates of the good data point
    i=coord.i(ind);
    j=coord.j(ind);
    
    % x-derivatives (We index by row, so if the depths are non zero, we know neighbors the are ind-1 and ind+1)
    if (H(i-1,j)~=0 && H(i+1,j)~=0)
        ind_W=ind-1;
        ind_E=ind+1;
        L1(ind,ind)=L1(ind,ind)-((H(i+1,j)+H(i,j))/2+(H(i,j)+H(i-1,j))/2)/(dx^2);
        L1(ind,ind_W)=L1(ind,ind_W)+(H(i,j)+H(i-1,j))/2/(dx^2);
        L1(ind,ind_E)=L1(ind,ind_E)+(H(i+1,j)+H(i,j))/2/(dx^2);
    elseif (H(i-1,j)==0 && H(i+1,j)~=0)
        ind_E=ind+1;
        L1(ind,ind)=L1(ind,ind)-(H(i+1,j)+H(i,j))/2/(dx^2);
        L1(ind,ind_E)=L1(ind,ind_E)+(H(i+1,j)+H(i,j))/2/(dx^2);
    elseif (H(i-1,j)~=0 && H(i+1,j)==0)
        ind_W=ind-1;
        L1(ind,ind)=L1(ind,ind)-(H(i,j)+H(i-1,j))/2/(dx^2);
        L1(ind,ind_W)=L1(ind,ind_W)+(H(i,j)+H(i-1,j))/2/(dx^2);
    end
    
    % y-derivatives (We have to find the north and south indicies)
    if (H(i,j-1)~=0 && H(i,j+1)~=0)
        ind_N=find(coord.i==i & coord.j==(j+1));
        ind_S=find(coord.i==i & coord.j==(j-1));
        L1(ind,ind)=L1(ind,ind)-((H(i,j+1)+H(i,j))/2+(H(i,j)+H(i,j-1))/2)/(dx^2);
        L1(ind,ind_S)=L1(ind,ind_S)+(H(i,j)+H(i,j-1))/2/(dx^2);
        L1(ind,ind_N)=L1(ind,ind_N)+(H(i,j+1)+H(i,j))/2/(dx^2);
    elseif (H(i,j-1)==0 && H(i,j+1)~=0)
        ind_N=find(coord.i==i & coord.j==(j+1));
        L1(ind,ind)=L1(ind,ind)-(H(i,j+1)+H(i,j))/2/(dx^2);
        L1(ind,ind_N)=L1(ind,ind_N)+(H(i,j+1)+H(i,j))/2/(dx^2);
    elseif (H(i,j-1)~=0 && H(i,j+1)==0)
        ind_S=find(coord.i==i & coord.j==(j-1));
        L1(ind,ind)=L1(ind,ind)-(H(i,j)+H(i,j-1))/2/(dx^2);
        L1(ind,ind_S)=L1(ind,ind_S)+(H(i,j)+H(i,j-1))/2/(dx^2);
    end
   
    PROGRESS_BAR(ind,1:N)
end

% Solve L1 eigenvalue problem  
[PHI lambda]=eigs(L1,Nm+1,'smallestabs');
lambda=diag(lambda);
lambda=-lambda(2:end);
PHI=PHI(:,2:end);

% Check L1 is self-adjoint
%PHI(:,2)'*(L1*PHI(:,1))./(PHI(:,1)'*(L1*PHI(:,1)))
%PHI(:,1)'*(L1*PHI(:,2))./(PHI(:,1)'*(L1*PHI(:,1)))

% Map and scale eigenvectors 
phi=zeros([Nx Ny Nm]);
mask=NaN([Nx Ny]);
mask(H>0)=1;

dA=dx^2;
A=nansum(nansum(mask.*dA));
for n=1:Nm
    tmp=zeros([Nx Ny]);
    tmp(water)=PHI(:,n);
    phi(:,:,n)=tmp;
    
    % Following Rao 1976
    tmp=lambda(n)/g*nansum(nansum(mask.*phi(:,:,n).^2.*dA))/A;
    phi(:,:,n)=phi(:,:,n)/sqrt(tmp);    
end

% Compute non-rotating periods
out.omega_i=sqrt(9.81*lambda);
out.T_i=2*pi./out.omega_i/(3600);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Define L2 matrix div(1/H grad()) with psi=0 on the boundary
water=find(H2>0); % Linear index to reconstruct modes
[coord.i,coord.j]=find(H2>0); % Row/column indicies to build matrix
N=length(coord.i);

% Set boundary to minimum depth for averages
H_tmp=H2;
H_tmp(H_tmp==0)=H_min; 

clear L1;
L2=sparse(N,N); 
for ind=1:N
    
    % The coordinates of the good data point
    i=coord.i(ind);
    j=coord.j(ind);
    
    L2(ind,ind)=L2(ind,ind)-2/(H_tmp(i+1,j)+H_tmp(i,j))/(dx^2) ...
                -2/(H_tmp(i,j)+H_tmp(i-1,j))/(dx^2);
    L2(ind,ind)=L2(ind,ind)-2/(H_tmp(i,j+1)+H_tmp(i,j))/(dx^2) ...
                -2/(H_tmp(i,j)+H_tmp(i,j-1))/(dx^2);
       
    if H2(i-1,j)~=0
        ind_W=ind-1;
        L2(ind,ind_W)=L2(ind,ind_W)+2/(H2(i,j)+H2(i-1,j))/(dx^2);
    end
    
    if H2(i+1,j)~=0
        ind_E=ind+1;
        L2(ind,ind_E)=L2(ind,ind_E)+2/(H2(i+1,j)+H2(i,j))/(dx^2);
    end
    
    if H2(i,j-1)~=0
        ind_S=find(coord.i==i & coord.j==(j-1));
        L2(ind,ind_S)=L2(ind,ind_S)+2/(H2(i,j)+H2(i,j-1))/(dx^2);
    end
    
    if H2(i,j+1)~=0
        ind_N=find(coord.i==i & coord.j==(j+1));
        L2(ind,ind_N)=L2(ind,ind_N)+2/(H2(i,j+1)+H2(i,j))/(dx^2);
    end
    
    PROGRESS_BAR(ind,1:N)
end

% Solve L2 eigenvalue problem  
[PSI mu]=eigs(L2,Nm,'smallestabs');
mu=-diag(mu);

% Check L2 is self-adjoint
%PSI(:,2)'*(L2*PSI(:,1))./(PSI(:,1)'*(L2*PSI(:,1)))
%PSI(:,1)'*(L2*PSI(:,2))./(PSI(:,1)'*(L2*PSI(:,1)))

% Map eigenvectors 
psi=zeros([Nx-1 Ny-1 Nm]);
mask=NaN([Nx-1 Ny-1]);
mask(H2>0)=1;

dA=dx^2;
A=nansum(nansum(mask.*dA));
for n=1:Nm
    tmp=zeros([Nx-1 Ny-1]);
    tmp(water)=PSI(:,n);
    psi(:,:,n)=tmp;
    
    % Following Rao 1976
    tmp=mu(n)/g*nansum(nansum(mask.*psi(:,:,n).^2.*dA))/A;
    psi(:,:,n)=psi(:,:,n)/sqrt(tmp);    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now find rotating gravity modes

% Compute transport for each mode on phi grid
M_phi=zeros([Nx Ny Nm]);
N_phi=zeros([Nx Ny Nm]);
M_psi=zeros([Nx Ny Nm]);
N_psi=zeros([Nx Ny Nm]);
for n=1:Nm
    for i=2:Nx-1
        for j=2:Ny-1
            % Phi volume flux
            if phi(i+1,j,n)~=0 && phi(i-1,j,n)~=0
                M_phi(i,j,n)=-H(i,j)*(phi(i+1,j,n)-phi(i-1,j,n))/(2*dx);
            end
            
            if phi(i,j+1,n)~=0 && phi(i,j-1,n)~=0
                N_phi(i,j,n)=-H(i,j)*(phi(i,j+1,n)-phi(i,j-1,n))/(2*dx);
            end
            
            % Psi volume flux
            M_psi(i,j,n)=-((psi(i-1,j,n)+psi(i,j,n))/2-(psi(i-1,j-1,n)+psi(i,j-1,n))/2)/dx;
            N_psi(i,j,n)=((psi(i,j,n)+psi(i,j-1,n))/2-(psi(i-1,j,n)+psi(i-1,j-1,n))/2)/dx;                 
        end
    end
end

% Compute coupling coefficients from integrals
mask=NaN([Nx Ny]);
mask(H>0)=1;
dA=dx^2;
area=nansum(nansum(mask.*dA));
for m=1:Nm
    for n=1:Nm
        A(m,n)=(nansum(nansum(mask.*f./(g*H).*(M_phi(:,:,m).*N_phi(:,:,n)-N_phi(:,:,m).*M_phi(:,:,n)).*dA)))/area;
        B(m,n)=(nansum(nansum(mask.*f./(g*H).*(M_phi(:,:,m).*N_psi(:,:,n)-N_phi(:,:,m).*M_psi(:,:,n)).*dA)))/area;
        D(m,n)=(nansum(nansum(mask.*f./(g*H).*(M_psi(:,:,m).*N_psi(:,:,n)-N_psi(:,:,m).*M_psi(:,:,n)).*dA)))/area;
    end
end
C=-B.';
nu=diag(sqrt(g*lambda));

% Now construct master matrix
Z=zeros([Nm Nm]);
a=[-A -B -nu; -C -D Z; nu Z Z];
%a=[-abs(A) -abs(B) -nu; -abs(C) -abs(D) Z; nu Z Z];

[PHI,sigma]=eig(ii*a);
sigma=diag(sigma);
sigma(sigma<1.1*f)=Inf;
[sigma,ind]=sort(sigma,1,'ascend');
PHI=PHI(:,ind);

% Assemble modes
P=PHI(1:Nm,:);
Q=PHI([1:Nm]+Nm,:);
R=PHI([1:Nm]+2*Nm,:);

R2=R(:,1:Nm).*repmat(sqrt(lambda/g),[1 Nm]);

Nm0=10;
eta=zeros([Nx Ny Nm0]);
M=zeros([Nx Ny Nm0]);
N=zeros([Nx Ny Nm0]);
for n=1:Nm0
    eta(:,:,n)=sum(repmat(permute(R2(:,n),[2 3 1]),[Nx Ny 1]).*phi,3);
    M(:,:,n)=sum(repmat(permute(P(:,n),[2 3 1]),[Nx Ny 1]).*M_phi+repmat(permute(Q(:,n),[2 3 1]),[Nx Ny 1]).*M_psi,3);
    N(:,:,n)=sum(repmat(permute(P(:,n),[2 3 1]),[Nx Ny 1]).*N_phi+repmat(permute(Q(:,n),[2 3 1]),[Nx Ny 1]).*N_psi,3);
end

% Normalize
mask=NaN([Nx Ny]);
mask(H>0)=1;
dA=dx^2;
area=nansum(nansum(mask.*dA));
amp=sqrt(real(squeeze(nansum(nansum((M.*conj(M)+N.*conj(N))./repmat(H,[1 1 Nm0])+g*eta.*conj(eta),1),2)*dA))/area);
for n=1:Nm0
    eta(:,:,n)=eta(:,:,n)/amp(n);
    M(:,:,n)=M(:,:,n)/amp(n);
    N(:,:,n)=N(:,:,n)/amp(n);
end

% Output
out.omega=sigma(1:Nm);
out.T=2*pi./out.omega/(3600);
out.phi=phi;
out.psi=psi;
out.eta=eta;
out.M=M;
out.N=N;
out.area=area;

return;
