istart=0;
idiag=0;

if (istart==1),
  clear all

  idiag=1;
  addpath /home/stephdut/matlab/m_map
  sizetyp=1;
% dirstr=['/nobackup1/stephdut/GUD_20151027/'];
% runstr=['run15_151_grazsame3']; nn=51;
% dirstr=['/nobackup1/stephdut/GUD_20170322/'];
% runstr=['run16_8']; nn=51;
% dirstr=['/nobackup1/stephdut/GUD_20180112/']; 
% runstr=['run17_29']; nn=51;
  dirstr=['/nobackup1/stephdut/GUD_20180501/'];
  runstr=['run18_30']; nn=52;

  time1=25920;


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


  gridname=['grid_ecco.nc']; 
  X=ncread(gridname,'X');
  Y=ncread(gridname,'Y');
  Z=ncread(gridname,'Z');
  drF=ncread(gridname,'drF');
  rA=ncread(gridname,'rA');
  HFacC=ncread(gridname,'HFacC');
  Zp1=ncread(gridname,'Zp1');
  Depth=ncread(gridname,'Depth');
 %
  nx=length(X); ny=length(Y); nz=length(Z);
  for k=1:nz,
    volc(k,:,:)=rA(:,:)*drF(k);
  end
  filand1=find(squeeze(HFacC(:,:,1))==0);
%%

  autoall=1; % 1=all, 2=only obligate autotrophs

  if (autoall==1), npl=35; else, npl=25; end
%%%
  for it=1:12,
    time=time1+240*it
    tstr=['00000',num2str(time)];
    for ityp=[1 3 4  5],
     if (ityp==1), ncid1=netcdf.open([run1str,'/3d.',tstr,'.nc']); str=['TRAC']; sstr=['phy']; ntyp=nn; end
     if (ityp==2), ncid1=netcdf.open([run1str,'/chl.',tstr,'.nc']); str=['TRAC']; sstr=['chl']; ntyp=npl; end
     if (ityp==3), ncid1=netcdf.open([run1str,'/PP.',tstr,'.nc']); str=['PP']; sstr=['pp']; ntyp=npl; end
     if (ityp==4), ncid1=netcdf.open([run1str,'/GR.',tstr,'.nc']); str=['GR']; sstr=['gr']; ntyp=nn; end

     if (ityp==5), ncid1=netcdf.open(['/home/stephdut/Diags_monod_radtrans9/NETCDF/mitgcm_ecco_theta_1992_2006.nc']); str=['THETA']; sstr=['theta']; ntyp=1; 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
%%%
      if (ityp==5),
        temp=squeeze(THETA(:,:,:,109:109+11));
      end % if
%%
     for ip=1:ntyp,
      if (ityp==1),
          eval([sstr,'(ip,:,:,:,it)=',str,num2str(ip+20),';']);
      end
      if (ityp==2),
          eval([sstr,'(ip,:,:,:,it)=',str,num2str(ip+20+nn),';']);
      end
      if (ityp==3 | ityp==4),
         if (ip<10), 
          eval([sstr,'(ip,:,:,:,it)=',str,'000',num2str(ip),';']);
         else
          eval([sstr,'(ip,:,:,:,it)=',str,'00',num2str(ip),';']);
         end
      end
     end % for ip
%%
     PPs(:,:,:,it)=PP;
     POCs(:,:,:,it)=TRAC12;
%
   end % for ityp
  end % for it
%%%%
 logvolexp=[-0.6900:.513:-0.6900+22*.513];
 logvol=(10*ones(size(logvolexp))).^logvolexp;
 logdm=2*((3/4*logvol/pi).^(1/3));
 if (nn==51),
    tsize=[1 2 3:4  5:9  5:9 5:15 7:16 7:22];
    ptype=[1 2 3 3 4 4 4 4 4 5 5  5 5 5 6 6 6 6 6  6 6 6 6 6 6 7 7 7  7 7 7 7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8];
 else
    tsize=[1 2 3:4  5:9  5:9 5:15 7:16 7:22 8];
    ptype=[1 2 3 3 4 4 4 4 4 5 5  5 5 5 6 6 6 6 6  6 6 6 6 6 6 7 7 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];
 end

QCARBON=[  2.5660936684839390E-012,  7.7891937212959131E-012,  2.3643540207836897E-011,  7.1768274555969635E-011,  2.1784746224399137E-010, 6.6126038419854616E-010,  2.0072085816665026E-009,  6.0927380296624019E-009,  1.8494070340848224E-008,  2.1784746224399137E-010,  6.6126038419854616E-010, 2.0072085816665026E-009,  6.0927380296624019E-009,  1.8494070340848224E-008,  2.1784746224399137E-010,  6.6126038419854616E-010,  2.0072085816665026E-009, 6.0927380296624019E-009,  1.8494070340848224E-008,  5.6137427230100969E-008,  1.7040114360624539E-007,  5.1724047885733581E-007,  1.5700464639297202E-006, 4.7657637011393349E-006,  1.4466134714414389E-005,  2.0072085816665038E-009,  6.0927380296624060E-009,  1.8494070340848238E-008,  5.6137427230101009E-008, 1.7040114360624550E-007,  5.1724047885733613E-007,  1.5700464639297210E-006,  4.7657637011393375E-006,  1.4466134714414398E-005,  4.3910916843308039E-005, 2.0072085816665038E-009,  6.0927380296624060E-009,  1.8494070340848238E-008,  5.6137427230101009E-008,  1.7040114360624550E-007,  5.1724047885733613E-007, 1.5700464639297210E-006,  4.7657637011393375E-006,  1.4466134714414398E-005,  4.3910916843308039E-005,  1.3328844615961171E-004,  4.0458754124946477E-004, 1.2280965323751224E-003,  3.7277991511405579E-003,  1.1315467591434890E-002,  3.4347292228348252E-002 6.092738029662405E-009];

%%
if (sizetyp==2),
  tsize(51)=8;
  QCARBON(51)=6.092738029662405E-009;
end % 
QCARBON=QCARBON';
%%%%
end % istart

PPtot=zeros(nx,ny,12);
GRtot=zeros(nx,ny,12);
clear tmp, tmp=squeeze(sum(gr,1));
clear tmp2, tmp2=squeeze(sum(pp,1));
%tmp2=pp; fi=find(pp~>0); %fi=find(isnan(pp)==1 | isinf(pp)==1);
%tmp2(fi)=0;
%tmp2=tmp2(1:35,:,:,:,:).*squeeze(phy(1:35,:,:,:,:));
%tmp2=squeeze(sum(tmp2,1));
for k=1:nz,
  PPtot(:,:,:)=PPtot(:,:,:)+squeeze(tmp2(:,:,k,:))*drF(k);
  GRtot(:,:,:)=GRtot(:,:,:)+squeeze(tmp(:,:,k,:))*drF(k);
end % for k

xtrans=-158;
GrX=180+(180+xtrans);
GrY=[22:1:50];
Grj=find(Y>=GrY(1)-0.5 & Y<=GrY(length(GrY))+0.5);
Gri=find(X>=GrX-0.5 & X<GrX+0.5);
GrYT=Y(Grj);

GrY2=[22:1:50];
Grj2=find(Y>=GrY2(1)-0.5 & Y<=GrY2(length(GrY2))+0.5);
GrYT2=Y(Grj2);

pltnt=0;
if (pltnt==1),
  kn1=1;
  figure
  typ=2; nnp=35; %nn;
  for it=1:4,
     if (it==1), mn1=2; cstr='c'; tstr=['Feb']; end
     if (it==2), mn1=4; cstr='r'; tstr=['Apr']; end
     if (it==3), mn1=6; cstr='k'; tstr=['Jun']; end
     if (it==4), mn1=8; cstr='m'; tstr=['Aug']; end
      gsst=squeeze(temp(Gri,Grj,1,mn1));
      clear fi, fi=find(gsst>18); t1=fi(length(fi)); lt1=GrYT(t1);
      clear fi, fi=find(gsst<13); if (isempty(fi)==0), t2=fi(1); lt2=GrYT(t2);
       else lt2=NaN; end
      for isize=1:3,
        clear fi,
        if (isize==1), fi=find(tsize(1:nnp)<5), scstr='g'; end 
        if (isize==2), fi=find(tsize(1:nnp)>4 & tsize(1:nnp)<8), scstr='b'; end        
        if (isize==3), fi=find(tsize(1:nnp)>9), scstr='m'; end
        bio=squeeze(sum(phy(fi,Gri,Grj,kn1,mn1),1))*12;
        if (typ==2), bio=bio./(squeeze(sum(phy(1:nnp,Gri,Grj,kn1,mn1),1))*12); end
        subplot(4,1,it)
        plot(GrYT,bio,[scstr,'-'],'LineWidth',2);
        hold on
      end % for
      if (typ==2), axis([22 49 0 1]); end
      c=axis,
      plot([lt1 lt1],[c(3) c(4)],'k--'); plot([lt2 lt2],[c(3) c(4)],'k--');
      axis([22 49 c(3) c(4)]);
      bolder
       if (typ==1), ylabel('biomass (mgC/m3)'); else, ylabel('fraction of biomass'); end
      text(23,c(3)+.8*(c(4)-c(3)),tstr);
  end % for it
  xlabel('latitude'); 
end % for ip

pltnt=0;
if (pltnt==1),
  kn1=1;
  figure
  typ=2; %1=absolute, 2=fraction
  for it=1:4,
     if (it==1), mn1=2; cstr='c'; tstr=['Feb']; end
     if (it==2), mn1=4; cstr='r'; tstr=['Apr']; end
     if (it==3), mn1=6; cstr='k'; tstr=['Jun']; end
     if (it==4), mn1=8; cstr='m'; tstr=['Aug']; end
      gsst=squeeze(temp(Gri,Grj,1,mn1));
      clear fi, fi=find(gsst>18); t1=fi(length(fi)); lt1=GrYT(t1);
      clear fi, fi=find(gsst<13); if (isempty(fi)==0), t2=fi(1); lt2=GrYT(t2);
       else lt2=NaN; end
      for isize=1:3,
        clear fi,
        if (isize==1), fi=find(tsize(1:35)<5), scstr='g'; end
        if (isize==2), fi=find(tsize(1:35)>4 & tsize(1:35)<8), scstr='b'; end
        if (isize==3), fi=find(tsize(1:35)>9), scstr='m'; end
        bio=squeeze(sum(pp(fi,Gri,Grj,kn1,mn1),1))*12*60*60*24;
        if (typ==2), bio=bio./(squeeze(sum(pp(:,Gri,Grj,kn1,mn1),1))*12*60*60*24); end
        subplot(4,1,it)
        plot(GrYT,bio,[scstr,'-'],'LineWidth',2);
        hold on
      end % for
      c=axis,
      plot([lt1 lt1],[c(3) c(4)],'k--'); plot([lt2 lt2],[c(3) c(4)],'k--');
      axis([22 49 c(3) c(4)]);
      bolder
      if (typ==1), ylabel('PP (mgC/m3/d)'); else, ylabel('fraction of PP'); end
      text(23,c(3)+.8*(c(4)-c(3)),tstr);
  end % for it
  xlabel('latitude');
end % for ip

pltnt=0;
if (pltnt==1),
  kn1=1;
  figure
  typ=1;
  for it=1:4,
     if (it==1), mn1=2; cstr='c'; tstr=['Feb']; end
     if (it==2), mn1=4; cstr='r'; tstr=['Apr']; end
     if (it==3), mn1=6; cstr='k'; tstr=['Jun']; end
     if (it==4), mn1=8; cstr='m'; tstr=['Aug']; end
      gsst=squeeze(temp(Gri,Grj,1,mn1));
      clear fi, fi=find(gsst>18); t1=fi(length(fi)); lt1=GrYT(t1);
      clear fi, fi=find(gsst<13); if (isempty(fi)==0), t2=fi(1); lt2=GrYT(t2);
       else lt2=NaN; end
      for isize=1:3,
        clear fi,
        if (isize==1), fi=find(tsize(1:25)<5), scstr='g'; end
        if (isize==2), fi=find(tsize(1:25)>4 & tsize(1:25)<8), scstr='b'; end
        if (isize==3), fi=find(tsize(1:25)>9), scstr='m'; end
        bio=squeeze(sum(phy(fi,Gri,Grj,kn1,mn1),1))*12*60*60*24;
        if (isize==1), fi=find(tsize(26:35)<5), scstr='g'; end
        if (isize==2), fi=find(tsize(26:35)>4 & tsize(26:35)<8), scstr='b'; end
        if (isize==3), fi=find(tsize(26:35)>9), scstr='m'; end
        if (isempty(fi)==0),bio2=squeeze(sum(phy(fi,Gri,Grj,kn1,mn1),1))*12*60*60*24;
          else, bio2=zeros(size(bio)); end


        subplot(4,1,it)
        plot(GrYT,bio,[scstr,'-'],'LineWidth',2);
        plot(GrYT,bio2,[scstr,'--'],'LineWidth',2);
        hold on
      end % for
      c=axis,
      plot([lt1 lt1],[c(3) c(4)],'k--'); plot([lt2 lt2],[c(3) c(4)],'k--');
      axis([22 49 c(3) c(4)]);
      bolder
      if (typ==1), ylabel('biomass (mgC/m3)'); else, ylabel('fraction of biomass'); end
      text(23,c(3)+.8*(c(4)-c(3)),tstr);
  end % for it
end % for ip




pltnt=0;
if (pltnt==1),
  figure
  typ=2;
  for it=1:4,
     if (it==1), mn1=2; cstr='c'; tstr=['Feb']; end
     if (it==2), mn1=4; cstr='r'; tstr=['Apr']; end
     if (it==3), mn1=6; cstr='k'; tstr=['Jun']; end
     if (it==4), mn1=8; cstr='m'; tstr=['Aug']; end
      gsst=squeeze(temp(Gri,Grj,1,mn1));
      clear fi, fi=find(gsst>18); t1=fi(length(fi)); lt1=GrYT(t1);
      clear fi, fi=find(gsst<13); t2=fi(1); lt2=GrYT(t2);
      if (typ==1), 
       gpp=squeeze(PPtot(Gri,Grj,mn1))*12*(60*60*24);
       ggr=squeeze(GRtot(Gri,Grj,mn1))*12*(60*60*24);
      end
      if (typ==2),
       gpp=squeeze(sum(pp(:,Gri,Grj,1,mn1),1))*12*(60*60*24);
       ggr=squeeze(sum(gr(:,Gri,Grj,1,mn1),1))*12*(60*60*24);
      end
     subplot(4,1,it)
     plot(GrYT,gpp,[cstr,'-'],'LineWidth',2);
     hold on
     plot(GrYT,-ggr,[cstr,':'],'LineWidth',2);
     plot(GrYT,gpp-ggr,[cstr,'--'],'LineWidth',2);
     c=axis,
     plot([lt1 lt1],[c(3) c(4)],'k--'); plot([lt2 lt2],[c(3) c(4)],'k--');
     bolder
     axis([22 49 c(3) c(4)]);
     if (typ==1), ylabel('mgC/m2/d'); end
     if (typ==2), ylabel('mgC/m3/d'); end
     text(23,c(3)+.8*(c(4)-c(3)),tstr);
     bolder
  end % for it
  xlabel('latitude')
end % for ip

pltnt=1;
if (pltnt==1),
  kn1=1;
  typ=4; % 1=PP full column; 2=CP full column; surface PP
  ng=2; nng=ng*2+1;
  figure
  for it=4:4,
%    if (it==1), mn1=2; cstr='c'; tstr=['Feb']; end
%    if (it==2), mn1=4; cstr='r'; tstr=['Apr']; end
%    if (it==3), mn1=6; cstr='k'; tstr=['Jun']; end
%    if (it==4), mn1=8; cstr='y'; tstr=['Aug']; end
     if (it==1), mn1=1; cstr='c'; tstr=['Jan']; end
     if (it==2), mn1=3; cstr='r'; tstr=['Mar']; end
     if (it==3), mn1=5; cstr='k'; tstr=['May']; end
     if (it==4), mn1=7; cstr='y'; tstr=['Jul']; end
%    if (it==1), mn1=9; cstr='c'; tstr=['Sep']; end
%    if (it==2), mn1=10; cstr='r'; tstr=['Oct']; end
%    if (it==3), mn1=11; cstr='k'; tstr=['Nov']; end
%    if (it==4), mn1=12; cstr='y'; tstr=['Dec']; end
     nGr=length(GrYT2);
     sbio(3:nGr+2)=GrYT2;
     pbio(3:nGr+2)=GrYT2;
      gsst=squeeze(temp(Gri,Grj2,1,mn1));
      clear fi, fi=find(gsst>18); t1=fi(length(fi)); lt1=GrYT(t1);
      clear fi, fi=find(gsst<13); if (isempty(fi)==0), t2=fi(1); lt2=GrYT(t2); else, lt2=NaN; end
      subplot(2,2,it)
      for isize=1:3,
        clear fi,
        sbio((isize-1)*nng+2:isize*nng+1,1)=isize;
        sbio((isize-1)*nng+2:isize*nng+1,2)=-(360-X(Gri-ng:Gri+ng));
        pbio((isize-1)*nng+2:isize*nng+1,1)=isize;
        pbio((isize-1)*nng+2:isize*nng+1,2)=-(360-X(Gri-ng:Gri+ng));
        if (isize==1), fi=find(tsize(1:35)<5), scstr='g'; end
        if (isize==2), fi=find(tsize(1:35)>4 & tsize(1:35)<8), scstr='b'; end
        if (isize==3), fi=find(tsize(1:35)>9), scstr='m'; end
        bio=squeeze(sum(phy(fi,[Gri-ng:Gri+ng],Grj2,kn1,mn1),1))*12;
        ppio=squeeze(sum(pp(fi,[Gri-ng:Gri+ng],Grj2,kn1,mn1),1))*12*(60*60*24);
        sbio((isize-1)*nng+2:isize*nng+1,3:nGr+2)=bio;
        pbio((isize-1)*nng+2:isize*nng+1,3:nGr+2)=ppio;
        if (typ==1)
         pptmp=squeeze(PPtot([Gri-ng:Gri+ng],Grj2,mn1))*12*(60*60*24);
         ystr=['PP (mgC/m2/d)'];
        end
        if (typ==2),
         pptmp=squeeze(PPtot([Gri-ng:Gri+ng],Grj2,mn1))*12*(60*60*24)- ...
               squeeze(GRtot([Gri-ng:Gri+ng],Grj2,mn1))*12*(60*60*24);
         ystr=['CP (mgC/m2/d)'];
        end % if
        if (typ==3)
         pptmp=squeeze(sum(pp(1:35,[Gri-ng:Gri+ng],Grj2,kn1,mn1),1))*12*(60*60*24);
         ystr=['PP (mgC/m3/d)'];
        end
        if (typ==4)
         pptmp=squeeze(sum(pp(1:35,[Gri-ng:Gri+ng],Grj2,kn1,mn1),1))*12*(60*60*24)- ...
               squeeze(sum(gr(1:35,[Gri-ng:Gri+ng],Grj2,kn1,mn1),1))*12*(60*60*24);
         ystr=['CP (mgC/m3/d)'];
        end
        if (isize==3), sbio((isize)*nng+2:(isize+1)*nng+1,3:nGr+2)=pptmp;
                       sbio((isize)*nng+2:(isize+1)*nng+1,1)=4;
                       sbio((isize)*nng+2:(isize+1)*nng+1,2)=-(360-X(Gri-ng:Gri+ng)); end
        plot(bio,pptmp,[scstr,'o'],'LineWidth',2);
        clear p s 
        [p s]=polyfit(bio,pptmp,1);
        [FF DELTA]=polyval(squeeze(p),bio,s);
        [r2 rmse]=rsquare(pptmp,FF);
        pval(it,isize,:)=p;
        r2val(it,isize)=r2;
        ffpl=polyval(p,[0:1:40],s); 
        hold on
        plot([0:1:40],ffpl,[scstr,'-']);
      end % for
      bolder
      bax=40;
      if (typ==1), axis([0 bax 0 600]), yp=50; end
      if (typ==2), axis([0 bax 0 500]), yp=50; end
      if (typ==3), axis([0 bax 0 40]),  yp=4; end
      if (typ==4), axis([0 bax 0 30]),  yp=3; end
      c=axis; 
%     plot([lt1 lt1],[c(3) c(4)],'k--'); plot([lt2 lt2],[c(3) c(4)],'k--');
%     axis([22 49 c(3) c(4)]);
      bolder
      if (it==1 | it==3), ylabel(ystr); end
      if (it==3 | it==4), xlabel('biomass (mgC/m3)'); end
      text(30,yp,[tstr]);  
  end % for it
end % for ip

