%tagliabue 
 load ../../data/obs/iron/tagliaua/taglia.asc
 D=taglia; %(1)month       (2)year (3)long (4)Lat  (5)depth (6)dFe (nM)

%geotrcaes
 geofil='../../data/obs/iron/GEOTRACES_IDP2017_v2/discrete_sample_data/netcdf/GEOTRACES_IDP2017_v2_Discrete_Sample_Data.nc';
 t   = ncread(geofil,'date_time')
 timeStr=datevec(datenum(0006,1,1)+t);
 yr = timeStr(:,1);
 mo= timeStr(:,2);
 lon  = ncread(geofil,'longitude');
 lat  = ncread(geofil,'latitude');
 dep  = ncread(geofil,'var2');
 dFe  = ncread(geofil,'var73'); 
 dFe2 = ncread(geofil,'var151');
 %expand data
  n=0;
  [I J]=size(dFe);
  for i=1:length(I)
   for j=1:length(J)
    n=n+1;
    YR(n)=yr(j); 
    LON(n)=lon(j);
    LAT(n)=lat(j);
    MO(n)=mo(j);
    DEP(n)=dep(j);
    DFE(n)=dFe(i,j);
    DFE2(n)=dFe2(i,j);
   end
  end
  dFe2=DFE2; dFe=DFE; dep=DEP; lat=LAT; lon=LON; mo=MO; yr=YR;;

%merge fish data with bottle data
 I=find(isfinite(dFe2)==1);
 for i=1:length(I)
  dFe(I(i))= meannan([dFe(I(i)) dFe2(I(i))]);
 end
 clear dFe2

%merge tagliabue with geotraces
 I=find(isfinite(dFe)==1);
 D=[D; mo(I) yr(I) lon(I) lat(I) dep(I) dFe(I)];

figure; hold on
 m_proj('Mercator')
 m_plot(D(:,3),D(:,4),'r.')
 m_coast
 m_grid
 title('Iron Observations')

figure;
 subplot(2,2,1)
  I=find(D(:,5)<100);
  hist(D(I,6))
  xlabel('dissolved iron concentration   [nM]')
  ylabel('absolute occurence')
  title('between 0 - 100m depth')
 subplot(2,2,2)
  I=find(D(:,5)>100 & D(:,5)<500);
  hist(D(I,6))
  xlabel('dissolved iron concentration   [nM]')
  ylabel('absolute occurence')
  title('between 100m - 500m depth')
 subplot(2,2,3)
  I=find(D(:,5)>500 & D(:,5)<2000);
  hist(D(I,6))
  xlabel('dissolved iron concentration   [nM]')
  ylabel('absolute occurence')
  title('between 500m - 2000m depth')
 subplot(2,2,4)
  I=find(D(:,5)>2000 & D(:,5)<5000);
  hist(D(I,6))
  xlabel('dissolved iron concentration   [nM]')
  ylabel('absolute occurence')
  title('between 2000m - 5000m depth')


%-------- publication ready
 I=find(D(:,5)<30 & D(:,4)<-30); in=[0:0.05:0.7];  col=jet(length(in));
figure; drawnow;
 m_proj('stereographic','lat',-90,'long',0,'radius',60)
 hold on

 
 for i=1:length(I)
  if D(I(i),6)>in(end)
   colI = length(in);
  else
   colI=interp1([in(1) in(end)],[1 length(in)],D(I(i),6),'linear');
  end
  [cc]=m_plot(D(I(i),3),D(I(i),4),'r.'); set(cc,'markersize',14,'color',col(round(colI),:));
 end


          amcolormap(jet,in); amcaxis(in);

          [cc]=amcolorbar('horiz',in);
          set(cc,'linewidth',2,'tickdir','out','fontweight','bold', ...
                  'xtick',in(1:2:end),'xticklabel',in(1:2:end),'position', ...
                  [0.130+0.15 0.1100+0.1 0.7750-0.3 0.0611/4], ...
                  'fontsize',10)

          m_grid('xtick',12,'tickdir','out','ytick',[-70:30:20],'linest','-', ...
                 'color','k','xticklabels',[],'yticklabels',[]);

          m_coast('patch',[0.7 0.7 0.7])
          legend('off')
  title('Iron Observations above 30m')
  brighten(0.2)
  print -dtiff ../ps/ironObsAll.tiff

%---seasonal variation
d=D(I,:); 
n=0;
for mo=1:12
 im=find(d(:,1)==mo);
 if isempty(im)~=1
  n=n+1;
  Fob(n)=mean(d(im,6)); 
  Fob_std(n)=std(d(im,6));
  T(n)=mo;
 end
end

figure
subplot(2,2,1)
set(gca,'nextplot','add','tickdir','out','linewidth',1.5, ...
        'fontsize',14,'box','on','xlim',[0.5 12.5],'ylim',[0 1.5 ])
%errplot_h(T,Fob,Fob_std,Fob_std)
%get model data
 Fmod=ncread('../../data/processed/feSeasoCy.nc','FES');
 Fmod_std=ncread('../../data/processed/feSeasoCy.nc','FESTD');

errplot_hb([1:12],Fmod,Fmod_std,Fmod_std)

errplot_h(T,Fob,Fob_std,Fob_std)

xlabel('month of year')
ylabel('dissolved Fe [nM]')
print -dpdf ../ps/ironSeonComp.pdf

%---vertical variation ----
figure; clear Fob, Fob_std
subplot(1,2,1)
set(gca,'nextplot','add','tickdir','out','linewidth',1.5, ... 
        'fontsize',14,'box','on','ydir','reverse','ylim',[0 5000],'xlim',[0 2])
%errplot_h(T,Fob,Fob_std,Fob_std)
%get model data
 Fmod=ncread('../../data/processed/feClima.nc','FES');
 Fmod_std=ncread('../../data/processed/feClima.nc','FESTD');
 z       =ncread('../../data/processed/feClima.nc','ST_OCEAN');
 zb=[(z(1:end-1)+z(2:end))/2]; zb=[0; zb; z(end)+125];

 n=0;
 for k=1:length(z)
  I=find(D(:,4)<-30 & D(:,5)>=zb(k) & D(:,5)<zb(k+1));
  if length(I)>0
   n=n+1;
   zob(n)=z(k);
   Fob(n)=mean(D(I,6));
   Fob_std(n)=std(D(I,6));
  end
 end

errplot_b(z,Fmod,Fmod_std,Fmod_std)
errplot(zob,Fob,Fob_std,Fob_std)
ylabel('depth   [m]')
xlabel('dissolved Fe [nM]')

print -dpdf ../ps/ironVertComp.pdf

