% uses mat files produced by read_car_nut
istart=1;
close all

if (istart==1),
 clear all
%
 logvolexp=[-1.41300:.513:-0.6900+25*.513];
 logvol=(10*ones(size(logvolexp))).^logvolexp;
 logdm=2*((3/4*logvol/pi).^(1/3));
%
 nphyt=31; nzoo=16; nbac=3; nplank=nphyt+nzoo+nbac;
%
 tsize=[2 3 4:5  6:10 6:10 7:15 8:15 7:22 1:3];
 ptype=[1 2 3 3 4 4 4 4 4 5 5  5 5 5 6 6 6 6 6  6 6 6 6  7 7 7  7 7 7 7 7  8 8 8 8 8 8 8 8 8 8 8 8 8 8  8 8 9 9 9];
%
 load grid_ecco.mat

 nx=length(X); ny=length(Y); nz=length(Z);
 for k=1:nz,
   volc(k,:,:)=rA(:,:)*drF(k);
 end
 ndep=10;
%
%% Gradients locations
 ix=180+(180-158); 
%
   GrY=[22:1:53]; GrY1=[22:1:54];

   Grj=find(Y>=GrY(1)-0.5 & Y<=GrY(length(GrY))+0.5);
   GrYT=Y(Grj);
%
 time=[1:12];
 for it=time,
   eval(['load /nfs/micklab001/stephdut/DAR3_MAT/run22_77_nobact_m',num2str(it),'_car_nut.mat']);
   for ip=1:nplank
     phy(it,ip,:,:)=squeeze(car(ip,ix,Grj,1:10));
     if (ip<nphyt+1), pps(it,ip,:,:)=squeeze(pp(ip,ix,Grj,1:10)); end
   end % for ip
   NO3(it,:,:)=squeeze(nutr(1,ix,Grj,1:10));
   DIN(it,:,:)=squeeze(sum(nutr(1:3,ix,Grj,1:10),1));
   PO4(it,:,:)=squeeze(nutr(4,ix,Grj,1:10));
   FET(it,:,:)=squeeze(nutr(5,ix,Grj,1:10));
   SIL(it,:,:)=squeeze(nutr(6,ix,Grj,1:10));
%
   T_DIN(it,:,:)=squeeze(sum(nutr_source(1:3,ix,Grj,1:10),1));
   T_PO4(it,:,:)=squeeze(nutr_source(4,ix,Grj,1:10));
   T_FET(it,:,:)=squeeze(nutr_source(5,ix,Grj,1:10));
   T_SIL(it,:,:)=squeeze(nutr_source(6,ix,Grj,1:10));
%
   SS_DIN(it,:,:)=squeeze(S_DIN(ix,Grj,1:10));
   SS_PO4(it,:,:)=squeeze(S_PO4(ix,Grj,1:10));
   SS_FET(it,:,:)=squeeze(S_Fe(ix,Grj,1:10));
   SS_SI(it,:,:)=squeeze(S_Si(ix,Grj,1:10));

 end % for it
%%%
 for ip=1:3,
    if (ip==1), ncid1=netcdf.open(['/home/stephdut/Diags_monod_radtrans9/NETCDF/mitgcm_ecco_phys_1992_2006.nc']); end
    if (ip==2), ncid1=netcdf.open(['/home/stephdut/Diags_monod_radtrans9/NETCDF/mitgcm_ecco_theta_1992_2006.nc']); end
    if (ip==3), ncid1=netcdf.open(['/home/stephdut/Diags_monod_radtrans9/NETCDF/mitgcm_ecco_salt_1992_2006.nc']); end
   [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);
   end % ip
   mld=squeeze(MLD(:,:,97:108));
   temp=squeeze(THETA(:,:,1:ndep,97:108));
   salt=squeeze(SALT(:,:,1:ndep,97:108));
%%
end % for istart


typ=1;

pltnt=1; tt=[1:12];
%% plot size classes along transect
if (pltnt==1)
 for it=tt,
  figure
  subplot(2,1,1),
  clear tmp
  if (typ==1), nsize=3; end
  if (typ==2), nsize=5; end
  if (typ==3), nsize=7; end
  if (typ==4), nsize=2; end
  if (typ==5), nsize=5; end

  for isize=1:nsize,
%% pro syn 1-9
   if (typ==1), 
    movfile='mov1_nb_stack_picosynpico_m';
    cax=20;
    fityp=find(ptype<8);
    if (isize==1), fi=find(tsize(fityp)==2); cstr='b'; end
    if (isize==2), fi=find(tsize(fityp)==3); cstr='c'; end
    if (isize==3), fi=find(tsize(fityp)>3 & tsize(fityp)<9); cstr='g'; end
   end % if
   if (typ==2)
%% pro syn 1-4, 4-20, 20-100
    cax=30;
    movfile='mov1_nb_stack_picosynsize_m';
    fityp=find(ptype<8);
    if (isize==1), fi=find(tsize(fityp)==2); cstr='b'; end
    if (isize==2), fi=find(tsize(fityp)==3); cstr='c'; end
    if (isize==3), fi=find(tsize(fityp)>3 & tsize(fityp)<7); cstr='g'; end
    if (isize==4), fi=find(tsize(fityp)>6 & tsize(fityp)<11); cstr='m'; end
    if (isize==5), fi=find(tsize(fityp)>10); cstr='y'; end
   end % if typ
%%%%
   if (typ==3),
%% pro syn 1-3, dino, cocco, diatom, diaz
    cax=30;
    movfile='mov1_stack_picosynpft_m';
    fityp=find(ptype<10);
    if (isize==1), fi=find(ptype==1); cstr='b'; iss=2; end
    if (isize==2), fi=find(ptype==2); cstr='c'; iss=2; end
    if (isize==3), fi=find(ptype==3); cstr='g'; iss=2; end
    if (isize==4), fi=find(ptype==7); cstr='m'; iss=2; end
    if (isize==5), fi=find(ptype==4); cstr='k'; iss=2; end
    if (isize==6), fi=find(ptype==6); cstr='y'; iss=2; end
    if (isize==7), fi=find(ptype==5); cstr='r'; iss=2; end
   end % if
%%
   if (typ==4),
% 1-9, diaz
    movfile='mov1_stack_diaz_m';
    cax=20;
    fityp=find(ptype<8);
    if (isize==1), fi=find(tsize(fityp)<9); cstr='g'; end
    if (isize==2), fi=find(ptype==5); cstr='r'; end
   end % if
%%%
%% pro syn het bact
   if (typ==5),
    movfile='mov1_stack_prosynhetbact_m';
    cax=20;
    fityp=find(ptype<10);
    if (isize==1), fi=find(ptype==1); cstr='b'; end
    if (isize==2), fi=find(ptype==2); cstr='c'; end
    if (isize==3), fi=find(ptype==9 & tsize==1); cstr='y'; end
    if (isize==4), fi=find(ptype==9 & tsize==2); cstr='y'; end
    if (isize==5), fi=find(ptype==9 & tsize==3); cstr='y'; end
   end % if
%%%
   if (isize==1),
       tmp=squeeze(sum(phy(it,fityp(fi),:,1),2))*12;
       tmpold=zeros(size(tmp));
   else
       tmpold=tmp;
       tmp=tmp+squeeze(sum(phy(it,fityp(fi),:,1),2))*12;
   end % if
%
   ny=length(tmp);
   yy=[Y(Grj)]; yy(ny+1:2*ny)=Y(flipud(Grj));
   xx=[tmpold]; xx(ny+1:2*ny)=flipud(tmp);
   fill(yy,xx,cstr)
   hold on
   bolder
  end % for isize
   axis([22 45 0 cax]);
%
   if (typ==1),
    text(23,cax-2,'1-9um','FontSize',12,'Color','g')
    text(23,cax-4,'Syn','Fontsize',12,'Color','c')
    text(23,cax-6,'Pro','FontSize',12,'Color','b')
   end % if
   if (typ==2),
    text(23,cax-2,'20-100um','FontSize',12,'Color','y')
    text(23,cax-5,'4-20um','FontSize',12,'Color','m')
    text(23,cax-8,'1-4um','FontSize',12,'Color','g')
    text(23,cax-11,'Syn','Fontsize',12,'Color','c')
    text(23,cax-14,'Pro','FontSize',12,'Color','b')
   end % if
   if (typ==3),
    text(23,cax-2,'diazotroph','FontSize',12,'Color','r')
    text(23,cax-5,'diatom','FontSize',12,'Color','y')
    text(23,cax-8,'cocco','FontSize',12,'Color','k')
    text(23,cax-11,'mixo dino','FontSize',12,'Color','m')
    text(23,cax-14,'1-4um','FontSize',12,'Color','g')
    text(23,cax-17,'Syn','Fontsize',12,'Color','c')
    text(23,cax-20,'Pro','FontSize',12,'Color','b')
   end % if
   if (typ==4),
    text(23,cax-4,'1-9um','FontSize',12,'Color','g')
    text(23,cax-2,'Diazotroph','FontSize',12,'Color','r')
   end % if
   if (typ==5),
    text(23,cax-2,'Hetrotrophic bacteria','FontSize',12,'Color','y')
    text(23,cax-4,'Syn','Fontsize',12,'Color','c')
    text(23,cax-6,'Pro','FontSize',12,'Color','b')
   end % if
%%
%  xlabel('latitude')
   ylabel('biomass (mg C/m3)');
   title(['MONTH ',num2str(it)]);
%%
  cax1=-0.01; cax2=0.3;
  subplot(2,1,2)
  for inut=1:3, 
   if (inut==3), tmp1=T_DIN; tmp2=SS_DIN; css='k'; ln=3; end
   if (inut==1), tmp1=T_PO4*16; tmp2=SS_PO4*10; css='m'; ln=2; end
   if (inut==2), tmp1=T_FET*1e4; tmp2=SS_FET*1e4; css='c'; ln=2; end
   tmp=squeeze(tmp1(it,:,1)+tmp2(it,:,1))*24*60*60;
   plot(Y(Grj),tmp,[css,'-'],'LineWidth',ln);
   hold on
   tmp=squeeze(tmp1(it,:,1))*24*60*60;
%  plot(Y(Grj),tmp,[csc,'--'],'LineWidth',2);
   tmp=squeeze(tmp2(it,:,1))*24*60*60;
%  plot(Y(Grj),tmp,[css,':'],'LineWidth',2);
  end
  axis([22 45 cax1 cax2]);
  bolder
  text(23,cax2-.02,'DIN','FontSize',12,'Color','k')
  text(23,cax2-.06,'16xPO4','Fontsize',12,'Color','m')
  text(23,cax2-.10,'1e4xFe','FontSize',12,'Color','c')

  xlabel('latitude'); ylabel('Supply Nutrients')
%
  if (it<10),
    eval(['print -dpng MOV/',movfile,'0',num2str(it),'.png']);
  else
    eval(['print -dpng MOV/',movfile,num2str(it),'.png']);
  end % if
 end % for it
end % for pltnt

