function plotWaterHeight(lambda, phi, cnm, snm, GM, R)

G   = 6.673e-11;
rho = 1025;   % density of sea water [kg/m**3]
maxDegree = size(cnm,1)-1;

%% read load love numbers
[n hn ln kn]  = textread('data/Load_Love2_CM.dat' , '%n%n%n%n' , 'headerlines',13);

%% compute factors for conversion of potential to equivalent water height
factor = zeros(1,maxDegree+1);
for n=2:1:maxDegree
  factor(n+1) = GM/(4*pi*G*rho*R^2)*(2*n+1)/(1+kn(n+1)); % for water height
  %factor(n+1) = GM/R; % for potential
end

%% compute values at grid
cosm   = cos([0:1:maxDegree]'*lambda);
sinm   = sin([0:1:maxDegree]'*lambda);
V = zeros(size(phi,2), size(lambda,2));
for i=1:size(phi,2)
  Pnm = legendreFunctions(pi/2-phi(i), maxDegree);
  V(i,:) = factor * ((cnm.*Pnm) * cosm + (snm.*Pnm) * sinm);
end
%% figure
contourf(lambda*180/pi, phi*180/pi, V, 70);
colorbar;
shading flat;
title('water height [m]','FontSize',16,'FontWeight','bold')
xlabel('Longitude','FontSize',12)
ylabel('Latitude','Fontsize',12)
axis tight
caxis([-0.5 0.5]);
