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));
%
%dirstr=['/nobackup1/stephdut/DAR7_20220808/DAR_EXAMPLES_1D/'];
%runstr=['run_6+4+0_RT_QUOTA_1D/']
%run1str=[dirstr,runstr,'/diags_20221021_0001'];
%dirstr=['/orcd/data/stephdut/001/stephdut/dar8/'];
 dirstr=['/home/stephdut/darwin3_20250122_vdarwin/OneD/'];
 runstr=['input_6+4+0_RT_QUOTA_1D/']
 run1str=[dirstr,runstr,'/diags_20250730_0001/'];
 runtyp=55; %49=31+16+3  %51=5+2+1 %52=5+2+2

%runstr=['run31_1_1D']
%runtyp=49; %49=31+16+3  %52=5+2+1
%runstr=['run_1D_421_bats']
%runtyp=54; %49=31+16+3  %52=5+2+1

nquota=1;
pquota=1;
fequota=1;

%%
%%
if (runtyp==49),
   nphyt=31; nzoo=16; nbac=3; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==51),
   nphyt=5; nzoo=2; nbac=1; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==52),
   nphyt=5; nzoo=2; nbac=2; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==53),
   nphyt=3; nzoo=2; nbac=1; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==54),
   nphyt=4; nzoo=2; nbac=1; nplank=nphyt+nzoo+nbac;
end % if
if (runtyp==55),
   nphyt=6; nzoo=4; 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,'/par.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='diag52b'; end
   if (runtyp==51), typstr='diag521'; end
   if (runtyp==52), typstr='diag522'; end
   if (runtyp==53), typstr='diag321'; end
   if (runtyp==54), typstr='diag421'; end
   if (runtyp==55), typstr='diag640'; 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_nit.data'],'%c','delimiter','\n','whitespace','');
       nplast=nplank;
       for ip=1:nplast
         str=mstr(6*(ip-1)+1:6*ip)';
         eval(['nit(',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
%
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=zeros(size(squeeze(PAR001)));
%for il=1:13,
%  if (il<10), eval(['tmp=PAR00',num2str(il),';']); end
%  if (il>9), eval(['tmp=PAR0',num2str(il),';']); end
%  totpar=totpar+squeeze(tmp);
%end % for il
totpar=squeeze(PAR);

end % for istart

 tax1=0; tax2=length(T);
%tax1=12; tax2=48;
%tax1=48; tax2=84;
dax2=-350;
cax1=-2; cax2=2.5;
cx=1;
%%
figure
 colormap(fake_parula())

for ip=1:6,
  subplot(6,1,ip),
  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]);
  title(['phytoplankton ',num2str(ip)]);
end % for


%%
figure
 colormap(fake_parula())

for ip=1:4,
  subplot(6,1,ip),
  tmp=squeeze(car(6+ip,:,:))*12;
  tmp=log10(tmp);
  pcolor(tim,-Z,tmp), 
  if (cx==1), caxis([cax1 cax2]); end
  shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
  title(['zooplankton ',num2str(ip)]);
end % for


figure
 colormap(hot)
 subplot(4,1,1);
 tmp=totpar;
  pcolor(tim,-Z,tmp), caxis([0 100]); shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
 title('PAR')

 subplot(4,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]);
 title('Chl');

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

 subplot(4,1,4)
 tmp=squeeze(TRAC12);
  pcolor(tim,-Z,tmp), shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
 title('POC');


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,tmp), shading flat, caxis([cax1 cax2]); colorbar
  axis([tax1 tax2 dax2 0]);
 title([tstr]);

end % for ip

figure
 colormap(fake_parula())
for ip=1:6,
  subplot(6,1,ip),
  tmp=squeeze(nit(ip,:,:));
  pcolor(tim,-Z,tmp),
  if (cx==1), caxis([0 0.1]); end
  shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
  title(['phyto ',num2str(ip),' N']);
end % for

figure
 colormap(hot)
for ip=1:6,
  subplot(6,1,ip),
  tmp=squeeze(car(ip,:,:))./squeeze(nit(ip,:,:));
  pcolor(tim,-Z,tmp),
  if (cx==1), caxis([4 16]); end
  shading flat, colorbar
  axis([tax1 tax2 dax2 0]);
  title(['phyto ',num2str(ip),' C:N']);
end % for

