lon_g = rdmds('../Grid/XG');
lat_g = rdmds('../Grid/YG');
lon_c = rdmds('../Grid/XC');
lat_c = rdmds('../Grid/YC');
height_c = rdmds('../Grid/RC');
height_f = rdmds('../Grid/RF');
%w=rdmds('W',NaN);

w_ave = rdmds('wVeltave',1752000);
u_ave = rdmds('uVeltave',1752000);
v_ave = rdmds('vVeltave',1752000);
%T_ave = rdmds('Ttave',13140000);
h=rdmds('../Depth/Depth');

%w_ave = (w_ave(:,:,:,4) + w_ave(:,:,:,5) +w_ave(:,:,:,6) + w_ave(:,:,:,7) +w_ave(:,:,:,8))./5;

%% Average vertical velocity
% Have to ignore nan-values in mean (don`t mean over z_values under the ground)
[m,n] = size(lat_c);
for i = 1:length(height_c)
    for j = 1:m
        for k = 1:n
            if -h(j,k)>=height_f(1,1,i)
                w_ave(j,k,i) = NaN;                   
            end
        end
    end
end

w_m = nanmean(w_ave,3);

% Only look at surface velocities
surf_u = squeeze(u_ave(:,:,1));
surf_v = squeeze(v_ave(:,:,1));

% Dont plot velocities over land
[r,c]= find(-h>=0);
for i = 1:length(r)
     surf_u(r(i),c(i)) = NaN;
     surf_v(r(i),c(i)) = NaN;
end

% To plot fewer arrows
new_lon=lon_g(1:6:end, 1:6:end);
new_lat=lat_g(1:6:end, 1:6:end);
n_surf_u = surf_u(1:6:end,1:6:end);
n_surf_v = surf_v(1:6:end,1:6:end);
nanre = ~isnan(n_surf_u) & ~isnan(n_surf_v);

V = [0, -500];

% newmap = jet;                    %starting map
% ncol = size(newmap,1);           %how big is it?
% zpos = 1 + floor(2/3 * ncol);    %2/3 of way through
% newmap(zpos,:) = [1 1 1];        %set that position to white
% colormap(newmap);                %activate it

%% Speed surface
% u_surf = squeeze(u_ave(:,:,1));
% v_surf = squeeze(v_ave(:,:,1));
% speed = sqrt(u_surf.^2 + v_surf.^2);
% 
% max_v = max(max(speed));
%% 
figure;  
pcolor(lon_g,lat_g,w_m);
caxis([-2*10^(-4) 1*10^(-4)])
colorbar
shading('interp')
freezeColors
hold on
contourf(lon_g,lat_g,-h,V,'color','k')
title('Surface velocities with pangean topography')

%nanre = ~isnan(surf_u) & ~isnan(surf_v);
%quiver(lon(nanre),lat(nanre),surf_u(nanre),surf_v(nanre))

quiver(new_lon(nanre),new_lat(nanre),n_surf_u(nanre),n_surf_v(nanre))

hold off

%new_u1=squeeze(u_ave1(1:4:end,1:3:end,1));
%new_v1=squeeze(v_ave1(1:4:end,1:3:end,1));

%figure; 
 
%pcolor(lon,lat,w_m1);
%caxis([-2*10^(-4) 0.1*10^(-4)])
%colorbar

%shading('interp')
%freezeColors
%hold on

%contourf(lon,lat,-h,V,'color','k')

%title('W average 300 year')

%quiver(new_lon,new_lat,new_u1,new_v1)
%hold off

%% make a sphere
% figure;
% [x,y,z] = sphere(31); 
% c = hadamard(32);
% colormap([1 0 1; 0 0 1])
% axis equal
% surf(x,y,z,c)


%% Change the contourcolour to pink:
% contourf(lon,lat,-h,[0 0], 'Color',[1,0.4,0.6]);

% Change background colour of figure:
% set(gcf,'Color',[1,0.4,0.6])
% also works with:
% greenbg([0 1 0]);


%% Vertical velocity different layers 
% layer = 13;
% str = 'W 100 year layer';
% figure;
% pcolor(lon,lat,w_ave2(:,:,layer));
% colorbar
% freezeColors
% shading('interp')
% hold on
% contourf(lon,lat,-h,24)
% title(sprintf('%s: %d', str, layer))
% hold off 

%% Vertical velocity different timesteps
% for i = 1:3
%     figure;
%     pcolor(lon,lat,w_m1(:,:,:,i));
%     colorbar
%     freezeColors
%     caxis([-1e-7 1e-7])
%     shading('interp')
%     hold on
%     contourf(lon,lat,-h,24)
% end

%% Temperature field
% 
% T_ave=squeeze(T_ave(:,:,1));
%  
% figure; 
% %set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1])
% pcolor(lon,lat,T_ave);
% colormap(jet)
% colorbar
% %freezeColors
% shading('interp')
% hold on
% contourf(lon,lat,-h,[0 0], 'Color', [1 1 1]);
% %colormap(gray)
% % Change the contourcolour to pink:
% % contourf(lon,lat,-h,[0 0], 'Color',[1,0.4,0.6]);
% title('T average 250 year')
% whitebg([1 1 1]);
%% Plotte færre piler: 
%  
% new_lon=lon(1:4:end, 1:3:end);
% new_lat=lat(1:4:end, 1:3:end);
% new_h=h(1:4:end, 1:3:end);
% new_u=squeeze(u_ave(1:4:end,1:3:end,1));
 %new_v=squeeze(v_ave(1:4:end,1:3:end,1));
 %new_u2=squeeze(u_ave(1:4:end,1:3:end,12));
 %new_v2=squeeze(v_ave(1:4:end,1:3:end,12));
% [m,n] = size(new_u);
% for i = 1:4:m
%     for j = 1:3:n
%         if (new_h(i,j) ~= 0)
%             quiver(new_lon(i,j),new_lat(i,j),new_u(i,j),new_v(i,j), 'Color', [0 0 1])
%             hold on 
%         end
%     end  
% end
%quiver(new_lon,new_lat,new_u,new_v, 'Color', 'k')
%hold on 
%quiver(new_lon,new_lat,new_u2,new_v2, 'Color', 'g')
%hold off
%% Horizontal velocities
% u_surf = squeeze(u_ave(:,:,1));
% v_surf = squeeze(v_ave(:,:,1));
% [m,n] = size(v_surf);
% for i = 1:m
%     for j = 1:n
%         if (h(i,j) ~= 0)
%             quiver(lon(i,j),lat(i,j),u_surf(i,j),v_surf(i,j), 'Color', [0 0 1])
%             hold on 
%         end
%     end  
% end
% hold off
% quiver(x,y,u_surf,v_surf)
%hold off

%% m_map:
% subplot(211);
% Slongs=[-100 0;-75 25;-5 45; 25 145;45 100;145 295;100 290];
% Slats= [  8 80;-80  8; 8 80;-80   8; 8  80;-80   0;  0  80];
% for l=1:7,
%  m_proj('sinusoidal','long',Slongs(l,:),'lat',Slats(l,:));
%  m_grid('fontsize',6,'xticklabels',[],'xtick',(-180:360),...
%         'ytick',(-90:90),'yticklabels',[],'linest','-','color',[.9 .9 .9]);
%  m_coast('patch','g');
% end;
% xlabel('Interrupted Sinusoidal Projection of World Oceans');
% % In order to see all the maps we must undo the axis limits set by m_grid calls:
% set(gca,'xlimmode','auto','ylimmode','auto');
% 
% subplot(212);
% Slongs=[-100 43;-75 20; 20 145;43 100;145 295;100 295];
% Slats= [  0  90;-90  0;-90   0; 0  90;-90   0;  0  90];
% for l=1:6,
%  m_proj('mollweide','long',Slongs(l,:),'lat',Slats(l,:));
%  m_grid('fontsize',6,'xticklabels',[],'xtick',(-180:30:360),...
%         'ytick',(-90:20:90),'yticklabels',[],'linest','-','color','k');
%  m_coast('patch',[.6 .6 .6]);
% end;
% xlabel('Interrupted Mollweide Projection of World Oceans');
% set(gca,'xlimmode','auto','ylimmode','auto');


