istart=1;
clear all
close all

 addpath /home/stephdut/matlab/m_map


for irun=1:4,
 if (irun==1)
%dutkiewicz et al 2015, BG (output from Diags_monod_radtrans9/plot_types_forzehr.m)
%%
  ncid1=netcdf.open(['dutkiewiczetal2015bg_depint.nc']);
  tric=1; unic=1; ddac=0; othc=0; conv=122; 
  grc=1;
 end % irun
 if (irun==3)
% follett  al 2018, ISME (output from Diags_monod_radtrans9/plot_types_dda.m)
%%
  ncid1=netcdf.open(['follettetal2018isme_depint.nc']);
  tric=1; unic=1; ddac=1; othc=0; conv=122;
  grc=1;
 end % irun
 if (irun==4)
% dutkiewicz et al 2020, GCB (output from Diags_redunt/plot_plankton_52, use run19_33_redunt_0
%%
  ncid1=netcdf.open(['dutkiewiczetal2020gcb_depint.nc']);
  tric=1; unic=1; ddac=0; othc=1; conv=1;
  grc=1;
 end % irun
if (irun==2)
% dutkiewicz et al 2019, NC (output from Diags_monod_radtrans9_igsm2/read_fordiaz.m)
%%
  ncid1=netcdf.open(['dutkiewiczetal2019nc_depint.nc']);
  tric=1; unic=1; ddac=0; othc=0; conv=122;
  grc=2;
 end % irun


%%%%

   [ndim,nvar,natt,unlim]=netcdf.inq(ncid1);

     for i=0:ndim-1,
      [dimname, dimlength]=netcdf.inqDim(ncid1,i);
      diml(i+1)=dimlength;
     end % for i

     for i=0:nvar-1;
      [varname, xtype,dimid,natt]=netcdf.inqVar(ncid1,i);
      if (length(dimid)==1),
        eval([varname,'=netcdf.getVar(ncid1,i,0,diml(dimid+1));']);
      end % if
      if (length(dimid)==2),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0],[diml(dimid(1)+1) diml(dimid(2)+1)] );']);
      end % if
      if (length(dimid)==3),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)] );']);
      end % if
      if (length(dimid)==4),
         eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)  diml(dimid(4)+1)] );']);
       end % if
        eval([varname,'=double(',varname,');']);
     end % for i
   netcdf.close(ncid1);
 if (grc==1); lat=latitude; long=longitude; land=LandMask; end
 X=longitude; Y=latitude;
 if (grc==2), % convert
   [XI YI]=meshgrid( X, Y);
   [XJ YJ]=meshgrid( long, lat);
   tmp=interp2(XI,YI,Diaz_unicell',XJ,YJ);
   clear Diaz_unicell;
   Diaz_unicell=tmp';
   tmp=interp2(XI,YI,Diaz_tricho',XJ,YJ);
   clear Diaz_tricho;
   Diaz_tricho=tmp';
   tmp=interp2(XI,YI,nonDiaz',XJ,YJ);
   clear nonDiaz;
   nonDiaz=tmp';
 end % id
%%
  nx=length(long); ny=length(lat);
%%
   totdiaz(irun,1:nx,1:ny)=0;
   if (tric==1),
     tricho(irun,:,:)=Diaz_tricho;
     totdiaz(irun,:,:)=squeeze(totdiaz(irun,:,:))+Diaz_tricho;
   else
     tricho(irun,1:nx,1:ny)=NaN;
   end
   if (unic==1),
     unicell(irun,:,:)=Diaz_unicell;
     totdiaz(irun,:,:)=squeeze(totdiaz(irun,:,:))+Diaz_unicell;
   else
     unicell(irun,1:nx,1:ny)=NaN;
   end
   if (ddac==1),
     dda(irun,:,:)=Diaz_dda;
     totdiaz(irun,:,:)=squeeze(totdiaz(irun,:,:))+Diaz_dda;
   else
     dda(irun,1:nx,1:ny)=NaN;
   end
   if (othc==1),
     other(irun,:,:)=Diaz_other;
     totdiaz(irun,:,:)=squeeze(totdiaz(irun,:,:))+Diaz_other;
   else
     other(irun,1:nx,1:ny)=NaN;
   end 
   nondiaz(irun,:,:)=nonDiaz;
%%
  figure
  cax1=-1; cax2=2.5;
  xwide=0.5; yhigh=0.25;
  colormap(fake_parula())
%
  for ityp=1:4,
   if (ityp==1), tmp=squeeze(tricho(irun,:,:)); end
   if (ityp==2), tmp=squeeze(unicell(irun,:,:)); end
   if (ityp==3), tmp=squeeze(dda(irun,:,:)); end
   if (ityp==4), tmp=squeeze(other(irun,:,:)); end
   tmp=log10(tmp);
%%
   ybottom=0.99-ityp*.25; xleft=0.2;
   h=axes('position',[(xleft+0.005) (ybottom) (xwide-0.01) (yhigh-0.01)]);
   projec=['mollweide'];
   m_proj(projec,'longitude',[0 360],'latitude',[-79 79]);
   m_pcolor(long,lat,tmp'), caxis([cax1 cax2]);
   shading flat,
   colorbar
   hp = gca ; hc = colorbar; axes(hc);
   bolder, set(gca,'Ytick',[-1 0  1 2] );
      set(gca,'YTickLabel',['0.1 ';'1.0 ';'10. ';'100.']);
   axes(hp)
   hold on
   m_coast('patch','k');
   m_grid('fontsize',6,'xticklabels',[],'xtick',[],...
       'ytick',[],'yticklabels',[],'linest','-','color',[0 0 0]);
   bolder
  end % for ityp

%%%
  figure(10)
 cax1=-1; cax2=2.5;
  xwide=0.5; yhide=0.25;
  colormap(fake_parula())
%
   tmp=squeeze(totdiaz(irun,:,:)); 
   tmp=log10(tmp);
%%
   ybottom=0.99-irun*.25; xleft=0.2;
   h=axes('position',[(xleft+0.005) (ybottom) (xwide-0.01) (yhigh-0.01)]);
   projec=['mollweide'];
   m_proj(projec,'longitude',[0 360],'latitude',[-79 79]);
   m_pcolor(long,lat,tmp'), caxis([cax1 cax2]);
   shading flat,
   colorbar
   hp = gca ; hc = colorbar; axes(hc);
   bolder, set(gca,'Ytick',[-1 0 1 2 ] );
      set(gca,'YTickLabel',['0.1 ';'1.0 ';'10. ';'100.']);
   axes(hp)
   hold on
   m_coast('patch','k');
   m_grid('fontsize',6,'xticklabels',[],'xtick',[],...
       'ytick',[],'yticklabels',[],'linest','-','color',[0 0 0]);
   bolder

end % for irun
writent=1;
if writent==1,
%
 ncid=netcdf.create(['foursimulations_depint.nc'],'NC_WRITE');
 dimnx=netcdf.defDim(ncid,'longitude',nx);
 dimny=netcdf.defDim(ncid,'latitude',ny);
 dimsim=netcdf.defDim(ncid,'simulations',irun);

 long_ID=netcdf.defVar(ncid,'longitude','double',[dimnx]);
 lat_ID=netcdf.defVar(ncid,'latitude','double',[dimny]);
 landmask_ID=netcdf.defVar(ncid,'LandMask','double',[dimnx,dimny]);
 Diazunicell_ID=netcdf.defVar(ncid,'Diaz_unicell','double',[dimsim,dimnx,dimny]);
 Diaztricho_ID=netcdf.defVar(ncid,'Diaz_tricho','double',[dimsim,dimnx,dimny]);
 Diazdda_ID=netcdf.defVar(ncid,'Diaz_dda','double',[dimsim,dimnx,dimny]);
 Diazother_ID=netcdf.defVar(ncid,'Diaz_other','double',[dimsim,dimnx,dimny]);
 nonDiaz_ID=netcdf.defVar(ncid,'nonDiaz','double',[dimsim,dimnx,dimny]);
 
 netcdf.endDef(ncid);
 
 netcdf.putVar(ncid,long_ID,long);
 netcdf.putVar(ncid,lat_ID,lat);

 netcdf.putVar(ncid,landmask_ID,land);
 netcdf.putVar(ncid,Diazunicell_ID,unicell);
 netcdf.putVar(ncid,Diaztricho_ID,tricho);
 netcdf.putVar(ncid,Diazdda_ID,dda);
 netcdf.putVar(ncid,Diazother_ID,other);
 netcdf.putVar(ncid,nonDiaz_ID,nondiaz);

 netcdf.close(ncid);
end % if

