% uses mat files produced by read_car_nut
istart=0;

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=[4 6 9];
 for it=time,
   eval(['load /nfs/micklab001/stephdut/DAR3_MAT/run22_77_m',num2str(it),'_car_nut.mat']);
   for ip=1:nphyt
     phy(it,ip,:,:)=squeeze(car(ip,ix,Grj,1:10));
     pps(it,ip,:,:)=squeeze(pp(ip,ix,Grj,1:10));
   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



%%%%%%%%
pltnt=0; tt=[4 6 9];
%% plot size classes along transect
if (pltnt==1)
 figure
 itt=0;
 for it=tt,
  itt=itt+1;
  subplot(3,1,itt),
  for isize=1:9,
   fityp=find(ptype<5 |ptype==6);
   fi=find(tsize(fityp)==1+isize); cstr='k-'; iss=2;
   if (isize==1), 
       tmp=squeeze(sum(phy(it,fityp(fi),:,1),2))*12;
   else
       tmp=tmp+squeeze(sum(phy(it,fityp(fi),:,1),2))*12;
   end % if
%
   plot(Y(Grj),tmp,cstr,'LineWidth',iss)
   hold on
   bolder
  end % for isize
   axis([22 45 0 10]);
   xlabel('latitude')
   ylabel('biomass (mg C/m3)');
 end % for it
end % for pltnt

pltnt=0; tt=[4 6 9];
%% plot size classes along transect
if (pltnt==1)
 figure
 itt=0;
 for it=tt,
  itt=itt+1;
  subplot(3,1,itt),
  clear tmp
  for isize=1:5,
   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
   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
   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 80]);
%  xlabel('latitude')
   ylabel('biomass (mg C/m3)');
 end % for it
end % for pltnt

pltnt=0;  tt=[4 6 9];
%% plot size classes along transect
if (pltnt==1)
 figure
 itt=0;
 for it=tt,
  itt=itt+1;
  subplot(3,1,itt),
  clear tmp
  for ip=1:7,
   if (ip==1), fi=find(ptype==1); cstr='b'; iss=2; end
   if (ip==2), fi=find(ptype==2); cstr='c'; iss=2; end
   if (ip==3), fi=find(ptype==3); cstr='g'; iss=2; end
   if (ip==4), fi=find(ptype==7); cstr='m'; iss=2; end
   if (ip==5), fi=find(ptype==4); cstr='w'; iss=2; end
   if (ip==6), fi=find(ptype==6); cstr='y'; iss=2; end
   if (ip==7), fi=find(ptype==5); cstr='r'; iss=2; end
   if (ip==1),
       tmp=squeeze(sum(phy(it,fi,:,1),2))*12;
       tmpold=zeros(size(tmp));
   else
       tmpold=tmp;
       tmp=tmp+squeeze(sum(phy(it,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 80]);
%  xlabel('latitude')
   ylabel('biomass (mg C/m3)');
 end % for it
end % for pltnt



%
%
pltnt=0; tt=[4 6 9];
%% plot size classes along transect
%% diaz/non diaz
%% functional types
if (pltnt==1)
 for it=tt,
  figure
  subplot(3,1,1), 
  for isize=1:3,
   if (isize==1), fi=find(tsize(1:nphyt)==2); cstr='g-'; iss=2; end
   if (isize==2), fi=find(tsize(1:nphyt)==3); cstr='y-'; iss=2; end
   if (isize==3), fi=find(tsize(1:nphyt)==4 | tsize(1:nphyt)==5); cstr='k-'; iss=2; end
   tmp=squeeze(sum(phy(it,fi,:,1),2))*12; 
%
   plot(Y(Grj),tmp,cstr,'LineWidth',iss)
   hold on
   bolder
   axis([22 45 0 10]);
   xlabel('latitude')
   ylabel('biomass (mg C/m3)');
  end % for isize

  subplot(3,1,2),
  for ip=1:2,
   if (ip==1), fi=find(ptype==5); cstr='c-'; iss=2; end
   if (ip==2), fi=find(ptype<5 | ptype==6 | ptype==7); cstr='k-'; iss=2; end
   tmp=squeeze(sum(phy(it,fi,:,1),2))*12;
%
   plot(Y(Grj),tmp,cstr,'LineWidth',iss)
   hold on
   bolder
   axis([22 45 0 10]);
   xlabel('latitude')
   ylabel('biomass (mg C/m3)');
  end % for ip
%%
  subplot(3,1,3),
  for ip=1:5,
   if (ip==1), fi=find(ptype<4); cstr='k-'; iss=2; end
   if (ip==2), fi=find(ptype==4); cstr='b-'; iss=2; end
   if (ip==3), fi=find(ptype==5); cstr='c-'; iss=2; end
   if (ip==4), fi=find(ptype==6); cstr='r-'; iss=2; end
   if (ip==5), fi=find(ptype==7); cstr='m-'; iss=2; end
   tmp=squeeze(sum(phy(it,fi,:,1),2))*12;
%
   plot(Y(Grj),tmp,cstr,'LineWidth',iss)
   hold on
   bolder
   axis([22 45 0 20]);
   xlabel('latitude')
   ylabel('biomass (mg C/m3)');
  end % for ip
%%
 end % for it
end % if


 
pltnt=1; tt=[4 6 9];
%%  plot nutrient sources
if (pltnt==1)
 titstr='Carbon Uptake'; cax1=0; cax2=0.5;
 figure
 itt=0;
 for it=tt,
  itt=itt+1;
% cmap = cmocean('balance'); colormap(cmap);
  colormap(fake_parula())
  tmp=squeeze(sum(pps(it,:,:,:),2))*24*60*60;
  subplot(3,1,itt), pcolor(Y(Grj),Z(1:10),tmp'); shading flat, caxis([cax1 cax2]);  colorbar
  hold on
% plot(Y(Grj),-mld(ix,Grj,it),'k-','LineWidth',2);
  axis([22 45 -300 0]);
  bolder
  ylabel('depth')
%
 end % for it
end % for pltnt



pltnt=0; tt=[4 6 9];
%%  plot nutrient sources
if (pltnt==1)
 for it=tt,
 figure
 for ip=1:3,
  if (ip==1), 
   tmp1=T_DIN; tmp2=SS_DIN; titstr='Supply of DIN'; cax1=-0.01; cax2=0.3;
  end
  if (ip==2), 
   tmp1=T_PO4; tmp2=SS_PO4; titstr='Supply of PO_4'; cax1=-0.001; cax2=0.02;
  end % if
  if (ip==3),
   tmp1=T_FET; tmp2=SS_FET; titstr='Supply of Fe'; cax1=-0.025*1e-4; cax2=0.2*1e-4;
  end
  subplot(3,1,ip)
  tmp=squeeze(tmp1(it,:,1)+tmp2(it,:,1))*24*60*60;
  plot(Y(Grj),tmp,'k-','LineWidth',2);
  hold on
  tmp=squeeze(tmp1(it,:,1))*24*60*60;
% plot(Y(Grj),tmp,'k--','LineWidth',2);
  tmp=squeeze(tmp2(it,:,1))*24*60*60;
% plot(Y(Grj),tmp,'k:','LineWidth',2);
  axis([22 45 cax1 cax2]);
  bolder
  xlabel('latitude'); ylabel(titstr)
 end % for ip
 end % for it
end % for pltnt

pltnt=0; tt=[4 6 9];
%%  plot nutrient sources
if (pltnt==1)
 for it=tt,
 figure
  cdiv=[18:4:59];
  colormap default, cmap=colormap; colormap(cmap(cdiv,:));

% cmap = cmocean('balance'); colormap(cmap);
%  colormap(fake_parula())
%
 for ip=1:3,
  if (ip==1),
   tmp1=squeeze(NO3(it,:,:)); titstr='NO3'; cax1=0; cax2=30;
  end
  if (ip==2),
   tmp1=squeeze(PO4(it,:,:)); titstr='PO_4'; cax1=0; cax2=2.5;
  end % if
  if (ip==3),
   tmp1=squeeze(FET(it,:,:)); titstr='Fe'; cax1=0; cax2=1e-3;
  end
%%
  subplot(3,1,ip), pcolor(Y(Grj),Z(1:10),tmp1'); shading flat, caxis([cax1 cax2]);  colorbar
  hold on
% plot(Y(Grj),-mld(ix,Grj,it),'k-','LineWidth',2);
  axis([22 45 -300 0]);
  bolder
  ylabel('depth')
  title(titstr);
 end % for ip
%%
 end % for it
end % for pltnt
