istart=1;
close all

if (istart==1),
clear all

% for run 31
 logvolexp=[-141300:.513:-0.6900+25*.513];
 logvol=(10*ones(size(logvolexp))).^logvolexp;
 logdm=2*((3/4*logvol/pi).^(1/3));
%
nquota=0;
pquota=0;
fequota=0;
%
 dirstr=['/orcd/data/stephdut/001/stephdut/HOME/darwin3_20240208/dar8/ONED/'];
%runstr=['run_2+1+0_1D/']; runtyp=56;
%runstr=['run_6+4+0_1D/']; runtyp=55;
%runstr=['run_6+4+0_RT_1D/']; runtyp=55;
 runstr=['run_6+4+0_RT_1D_NAT/']; runtyp=55;
%runstr=['run_6+4+0_RT_QUOTA_1D/']; runtyp=55; nquota=1; pquota=1; fequota=1;
%runstr=['run_31+16+3_RT_1D/']; runtyp=49;

 run1str=[dirstr,runstr,'/diags_20240329_0002'];


%%
%%
if (runtyp==49),
   nphyt=31; nzoo=16; nbac=3; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==50),
   nphyt=31; nzoo=32; nbac=3; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==55),
   nphyt=6; nzoo=4; nbac=0; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==56),
   nphyt=2; nzoo=1; nbac=0; nplank=nphyt+nzoo+nbac;
end % if




 numquota=nquota+pquota+fequota;
 ndiag=4+nquota+pquota+fequota;
 for idiag=1:ndiag,
   if (idiag==1), ncid1=netcdf.open([run1str,'/biomass_c.0000000000.t001.nc']); end
   if (idiag==2), ncid1=netcdf.open([run1str,'/chl.0000000000.t001.nc']); end
   if (idiag==3), ncid1=netcdf.open([run1str,'/General.0000000000.t001.nc']); end
   if (idiag==4), ncid1=netcdf.open([run1str,'/General_2D.0000000000.t001.nc']); end
   if (idiag==5),  ncid1=netcdf.open([run1str,'/biomass_n.0000000000.t001.nc']); end
   if (idiag==6),  ncid1=netcdf.open([run1str,'/biomass_p.0000000000.t001.nc']); end
   if (idiag==7),  ncid1=netcdf.open([run1str,'/biomass_fe.0000000000.t001.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
%%%
 end % for idiag

   if (runtyp==49), typstr='diag31163'; end
   if (runtyp==50), typstr='diag31323'; end
   if (runtyp==55), typstr='diag640'; end
   if (runtyp==56), typstr='diag210'; end
       clear mstr
       mstr=textread([typstr,'_list_rad_car.data'],'%c','delimiter','\n','whitespace','');
       nplast=nplank; 
       for ip=1:nplast
         str=mstr(6*(ip-1)+1:6*ip)';
         eval(['car(',num2str(ip),',:,:,:)=',str,';']);
       end % for ip
% 
       clear mstr
         mstr=textread([typstr,'_list_rad_chl_',num2str(numquota),'quota.data'],'%c','delimiter','\n','whitespace','');
         nplast=nphyt;
         for ip=1:nplast
          str=mstr(6*(ip-1)+1:6*ip)';
          eval(['chl(',num2str(ip),',:,:,:)=',str,';']);
         end % for ip

     if (nquota==1),
       clear mstr
         mstr=textread([typstr,'_list_rad_nit.data'],'%c','delimiter','\n','whitespace','');
         nplast=nplank;
         for ip=1:nplast
          str=mstr(6*(ip-1)+1:6*ip)';
          eval(['nbio(',num2str(ip),',:,:,:)=',str,';']);
         end % for ip
      end % if nquota

     if (pquota==1),
       clear mstr
         mstr=textread([typstr,'_list_rad_phos.data'],'%c','delimiter','\n','whitespace','');
         nplast=nplank;
         for ip=1:nplast
          str=mstr(6*(ip-1)+1:6*ip)';
          eval(['pbio(',num2str(ip),',:,:,:)=',str,';']);
         end % for ip
      end % if pquota

     if (fequota==1),
       clear mstr
         mstr=textread([typstr,'_list_rad_fet.data'],'%c','delimiter','\n','whitespace','');
         nplast=nplank;
         for ip=1:nplast
          str=mstr(6*(ip-1)+1:6*ip)';
          eval(['febio(',num2str(ip),',:,:,:)=',str,';']);
         end % for ip
      end % if fequota



%
tim=[1:1:length(T)];
delZ=[10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, 10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, 31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, 93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, 139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, 341.50,364.50,387.50,410.50,433.50,456.50];
 nz=length(delZ);
Z(1)=delZ(1)/2;
for k=2:nz
 Z(k)=Z(k-1)+delZ(k-1)/2+delZ(k)/2;
end 

totpar=squeeze(PAR);

end % for istart

%tax1=0; tax2=length(T); nint=48;
 tax1=0; tax2=4*12; nint=12;
%tax1=48; tax2=96; nint=12;

dax2=-350;
cax1=-2; cax2=2.5;
cx=1;
%%
% plot plankton %%%
if (runtyp==56),
   npanel=6; nfig=2;  ncol=1; np1=[1 3]; np2=[2 3];
   titlestr=['phyto';'zoo  '];
end
if (runtyp==55), 
   npanel=6; nfig=2; ncol=1; np1=[1 7]; np2=[6 10];  
   titlestr=['phyto';'zoo  '];
end 
if (runtyp==49), 
   npanel=5; nfig=7; ncol=3; 
%  np1=[1 5 10 15 20 24 29 32 37 42 47 48]; 
%  np2=[4 9 14 19 23 28 31 36 41 46 47 50]; 
%  titlestr=['pico  ';'cocco ';'diazo ';'diatom';'diatom';'mixo  ';'mixo  ';'zoo1-5';'z6-10 ';'z11-15';'zoo16 ';'bact  ']; 
   np1=[1 5 10 15 24 32 48 ];
   np2=[4 9 14 23 31 46 50 ];
   titlestr=['pico  ';'cocco ';'diazo ';'diatom';'mixo  ';'zoo   ';'bact  '];

end
if (runtyp==50),
   npanel=5; nfig=8; ncol=3;
%  np1=[1 5 10 15 20 24 29 32 37 42 47 48 52 57 62 64];
%  np2=[4 9 14 19 23 28 31 36 41 46 47 51 56 61 62 66];
%  titlestr=['pico  ';'cocco ';'diazo ';'diatom';'diatom';'mixo  ';'mixo  ';'zoo1-5';'z6-10 ';'z11-15';'Zoo16 ';'Zoo1-5';'Z6-10 ';'Z11-15';'Zoo16 ';'bact  '];
   np1=[1 5 10 15 24 32 48 64];
   np2=[4 9 14 23 31 46 61 66];
   titlestr=['pico  ';'cocco ';'diazo ';'diatom';'mixo  ';'zoo1  ';'Zoo2  ';'bact  '];
end

for ifig=1:nfig,
 figure
 colormap(fake_parula())

 ip1=0;
 for ip=np1(ifig):np2(ifig),
  ip1=ip1+1;
  subplot(npanel,ncol,ip1),
  tmp=squeeze(car(ip,:,:))*12;
  tmp=log10(tmp);
  pcolor(tim,-Z,tmp), 
  if (cx==1), caxis([cax1 cax2]); end
  shading flat, %colorbar
  axis([tax1 tax2 dax2 0]);
  set(gca,'Xtick',[0:nint:tax2]),
  if (ip1==1), title(titlestr(ifig,:)); end
 end % for
%xlabel('months'); ylabel('depth');
end % for ifig

if (nquota==1 | pquota==1 | fequota==1),
 whichquota=3; %1=nitrogen, 2=phos, 3=fe;
 if (whichquota==1), tmpq=nbio; cax1=4; cax2=20; end
 if (whichquota==2), tmpq=pbio; cax1=40; cax2=200; end
 if (whichquota==3), tmpq=febio; cax1=5e6; cax2=.8e7; end
  for ifig=1:nfig,
 figure
 colormap('spring')

 ip1=0;
 for ip=np1(ifig):np2(ifig),
  ip1=ip1+1;
  subplot(npanel,ncol,ip1),
  tmp1=squeeze(tmpq(ip,:,:));
  tmp2=squeeze(car(ip,:,:));
  fi=find(tmp2<1e-4);
  tmp=tmp2./tmp1;
  tmp(fi)=NaN;
  pcolor(tim,-Z,tmp),
  if (cx==1), caxis([cax1 cax2]); end
  shading flat,  colorbar
  axis([tax1 tax2 dax2 0]);
  set(gca,'Xtick',[0:nint:tax2]),
  if (ip1==1), title(titlestr(ifig,:)); end
 end % for
%xlabel('months'); ylabel('depth');
 end % for ifig
end % if

%%
figure
 colormap(hot)
 subplot(5,1,1);
 tmp=totpar;
  pcolor(tim,-Z,tmp), caxis([0 100]); shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
   set(gca,'Xtick',[0:nint:tax2]);
 title('PAR')

 subplot(5,1,2)
 tmp=squeeze(sum(chl,1));
 tmp=log10(tmp);
  pcolor(tim,-Z,tmp), caxis([-2 0.5]);  shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
   set(gca,'Xtick',[0:nint:tax2]),
 title('Chl');

 subplot(5,1,3)
 tmp=squeeze(TRAC20);
  pcolor(tim,-Z,tmp), shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
   set(gca,'Xtick',[0:nint:tax2]),
 title('CDOM');

 subplot(5,1,4)
 tmp=squeeze(TRAC08);
  pcolor(tim,-Z,tmp), shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
   set(gca,'Xtick',[0:nint:tax2]),
 title('DOC');
  set(gca,'Xtick',[0:nint:tax2]),

 subplot(5,1,5)
 tmp=squeeze(TRAC12);
  pcolor(tim,-Z,tmp), shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
 title('POC');
  set(gca,'Xtick',[0:nint:tax2]),
  xlabel('months'); ylabel('depth'); 


figure
colormap(copper)
for ip=1:4,
 if (ip==1), tmp=squeeze(TRAC02); tstr='NO3';
         cax1=0; cax2=5; end
 if (ip==2), tmp=squeeze(TRAC05); tstr='PO4'; 
         cax1=0; cax2=0.5; end
 if (ip==3), tmp=squeeze(TRAC06); tstr='Fe'; 
         cax1=0; cax2=1e-5; end
 if (ip==4), tmp=squeeze(TRAC07);  tstr='SiO2'; 
         cax1=0; cax2=2; end 
 subplot(4,1,ip)
  pcolor(tim,-Z(1:21),squeeze(tmp(1:21,:))), shading flat, %caxis([cax1 cax2]); 
  colorbar
  axis([tax1 tax2 dax2 0]);
   set(gca,'Xtick',[0:nint:tax2]),
 title([tstr]);
 if (ip==4), xlabel('months'); ylabel('depth'); end
end % for ip

pltnt=1;
if (pltnt==1),
figure
cmap=cmocean('thermal'); colormap(cmap)
for ip=1:2,
 if (ip==1), tmp=squeeze(THETA); tstr='temperature';
         cax1=8; cax2=15; end
 if (ip==2), tmp=squeeze(KPPdiffS); tstr='KPP diffusion';
         cax1=-0.05; cax2=0.05; end
 subplot(3,1,ip)
  pcolor(tim,-Z(1:21),squeeze(tmp(1:21,:))), shading flat, %caxis([cax1 cax2]); 
  colorbar
  hold on
  plot(tim,-squeeze(MXLDEPTH),'k-')
  axis([tax1 tax2 dax2 0]);
  set(gca,'Xtick',[0:nint:tax2]),
  if (ip==2), xlabel('month'); end
 title([tstr]);

end % for ip
end % pltnt

