function sphcolorcontdoubledouble(field1a,field1b,field1c,field2a,field2b,field2c,lat,lon,num)

%compares two color and contour spherical plots
%is written with the Atlantic domain in mind

%field1a-subplot 1 color
%field1b- subplot 1 contour





%set default figure n




%plots spherical data on a robinson projection with continents

if max(lon)<=4
    error('input longitude should be in degrees')
end

if max(abs(lat))>90
    error('input order is wrong - input latitude vector before longitude vector')
end

field1a=squeeze(field1a);
[NN,N]=size(field1a);



if length(lat)<NN | length(lat)>NN
    error('field1a size and latitude input must be same size')
end

if length(lon)<N | length(lon)>N
    error(' field1a size and longitude input must be same size')
end

field1b=squeeze(field1b);
[NN,N]=size(field1b);

if length(lat)<NN | length(lat)>NN
    error('field 1b size and latitude input must be same size')
end

if length(lon)<N | length(lon)>N
    error(' field 1b size and longitude input must be same size')
end

field2a=squeeze(field2a);
[NN,N]=size(field2a);



if length(lat)<NN | length(lat)>NN
    error('field 2a size and latitude input must be same size')
end

if length(lon)<N | length(lon)>N
    error(' field 2a size and longitude input must be same size')
end

field2b=squeeze(field2b);
[NN,N]=size(field2b);

if length(lat)<NN | length(lat)>NN
    error('field 2b size and latitude input must be same size')
end

if length(lon)<N | length(lon)>N
    error(' field 2b size and longitude input must be same size')
end

%convert longitudes to be continous
for n=1:N
    if lon(n) >= 180
        lon(n)=lon(n)-360;
    else
        lon(n)=lon(n);
    end
end


LAT=lat;
load coast
figure(num); clf; subplot(2,2,1)


axesm('robinson','MapLatLimit',[min(lat) max(lat)],'MapLonLimit',[min(lon) max(lon)], 'Origin', 0, 'FLatLimit',[min(LAT)-5 max(LAT)+5], 'FLonLimit',[min(lon) max(lon)],'ScaleFactor',1.5)

framem; plotm(lat,long, 'k');
Lat=lat;
lat=LAT; clear LAT



[LON,LAT]=meshgrid(lon,lat);
hold on;
cmap = [ ...
    1.0000    1.0000    1.0000
    0.6667    1.0000    0.6667
    0.3333    1.0000    0.3333
         0    0.9192         0
    1.0000    1.0000         0
    1.0000    0.5000         0
    0.8000         0         0
];
colormap(cmap);
caxis([0 max(max(abs([field1a field2a])))])
pcolorm(LAT, LON, field1a)
colorbar;
hold on
[c,h]=contourm(LAT, LON, field1b, 'm');
ht = clabelm(c,h); set(ht,'color','m')
hidem(gca)

hold on;framem; plotm(Lat,long, 'k');
tightmap

%legend('LGM')
title('LGM')

figure(num); hold on; subplot(2,2,3)


axesm('robinson','MapLatLimit',[min(lat) max(lat)],'MapLonLimit',[min(lon) max(lon)], 'Origin', 0, 'FLatLimit',[min(lat)-5 max(lat)+5], 'FLonLimit',[min(lon) max(lon)],'ScaleFactor',1.5)

framem; plotm(Lat,long, 'k');



[LON,LAT]=meshgrid(lon,lat);
hold on;
colormap(cmap);
caxis([0 max(max(abs([field1a field2a])))])
pcolorm(LAT, LON, field2a)
%colorbar;
hold on
[c,h]=contourm(LAT, LON, field2b, 'm');
ht = clabelm(c,h); set(ht,'color','m')
hidem(gca)

hold on;framem; plotm(Lat,long, 'k');
tightmap

%legend('MODERN');

colorbar

LAT=lat;
load coast
figure(num); hold on; subplot(2,2,2)


axesm('robinson','MapLatLimit',[min(lat) max(lat)],'MapLonLimit',[min(lon) max(lon)], 'Origin', 0, 'FLatLimit',[min(LAT)-5 max(LAT)+5], 'FLonLimit',[min(lon) max(lon)],'ScaleFactor',1.5)

framem; plotm(lat,long, 'k');
Lat=lat;
lat=LAT; clear LAT



[LON,LAT]=meshgrid(lon,lat);
hold on;
colormap(cmap);
caxis([0 max(max(abs([field1c field2c])))])
pcolorm(LAT, LON, field1c)
colorbar;
hold on
[c,h]=contourm(LAT, LON, field1b, 'm');
ht = clabelm(c,h); set(ht,'color','m')
hidem(gca)

hold on;framem; plotm(Lat,long, 'k');
tightmap

%legend('LGM')
title('LGM')

figure(num); hold on; subplot(2,2,4)


axesm('robinson','MapLatLimit',[min(lat) max(lat)],'MapLonLimit',[min(lon) max(lon)], 'Origin', 0, 'FLatLimit',[min(lat)-5 max(lat)+5], 'FLonLimit',[min(lon) max(lon)],'ScaleFactor',1.5)

framem; plotm(Lat,long, 'k');



[LON,LAT]=meshgrid(lon,lat);
hold on;
colormap(cmap);
caxis([0 max(max(abs([field1c field2c])))])
pcolorm(LAT, LON, field2c)
%colorbar;
hold on
[c,h]=contourm(LAT, LON, field2b, 'm');
ht = clabelm(c,h); set(ht,'color','m')
hidem(gca)

hold on;framem; plotm(Lat,long, 'k');
tightmap

%legend('MODERN');

colorbar

% figure(num); hold on; subplot(3,1,3)
% 
% 
% axesm('robinson','MapLatLimit',[min(lat) max(lat)],'MapLonLimit',[min(lon) max(lon)], 'Origin', 0, 'FLatLimit',[min(lat)-5 max(lat)+5], 'FLonLimit',[min(lon) max(lon)],'ScaleFactor',1.5)
% 
% framem; plotm(Lat,long, 'k');
% 
% 
% 
% [LON,LAT]=meshgrid(lon,lat);
% hold on;
% color_trans; colormap(cmap);
% caxis([-max(max(abs([field1a field2a]))) max(max(abs([field1a field2a])))])
% pcolorm(LAT, LON, field1a-field2a)
% %colorbar;
% hold on
% [c,h]=contourm(LAT, LON, field2b, 'r');
% ht = clabelm(c,h); set(ht,'color','r')
% hidem(gca)
% 
% hold on;framem; plotm(Lat,long, 'k');
% tightmap
% 
% legend('LGM-MODERN')



