clear all
close all

nplank=2; 
%mum=[.6 .9]; %1/d
%kr=[0.005 0.01]; %mmolN/m3
 mum=[3 2];
 kr=[0.05 0.07];
mp=[.1 .1]; %1/d
grm=[5 5 ]; %m3/(mmolN d)
ass=[.6 .6]; % unitless
mz=[.1 .1];  % 1/d
%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%
mucrit=(mum(2)*kr(1))/kr(2) + (mp(1)*(kr(2)-kr(1)))/(kr(2));
%%%%%%
Sr=[1e-4:5e-4:1e-2 1.1e-2:1e-3:2e-1]; % 1e-3+1e-3:1e-3:.5e-2];  %mmolN/(m3 d)
%%%%%
AA=(mp(1)/grm(1)) * (mum(1)*kr(2) - mum(2)*kr(1) - mp(1)*(kr(2)-kr(1)))/ ...
              (kr(1)*mum(2)+mp(1)*(kr(2)-kr(1)));
CC=AA; %???
srcrit(1)=(mp(1)*mz(1))/(ass(1)*grm(1));
srcrit(2)=srcrit(1)+(mz(1)*AA)/ass(1);
srcrit(3)=srcrit(1)+(mp(2)*mz(2))/(ass(2)*grm(2))+(mz(1)*CC)/ass(1);
srmu1crit(1)=100;
srmu1crit(2)=(mum(1)*mz(1))/(ass(1)*grm(1));
srmu1crit(3)=100;
srmu1crit(4)=100;
z1crit(1)=0;
z1crit(2)=(mum(1)+mp(1))/grm(1);
%%%%%
% REGION 1
clear fi, fi=find(Sr<srcrit(1));
Rs(fi)=(kr(1)*mp(1))./(mum(1)-mp(1));
Bs(fi,1)=Sr(fi)/mp(1);
Zs(fi,1)=0;
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=Rs(fi);

% REGION 2
clear fi, fi=find(Sr>srcrit(1) & Sr<=srcrit(2) & Sr<srmu1crit(2));
Rs(fi)=(kr(1)*ass(1)*grm(1)*Sr(fi))./(mum(1)*mz(1)-ass(1)*grm(1)*Sr(fi));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Zs(fi,1)=(Sr(fi)*ass(1))/mz(1) - mp(1)/grm(1);
Bs(fi,2)=0;
Zs(fi,2)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));

clear f1i, f1i=find(Zs(:,1)>z1crit(2));

% REGION 3
clear fi, fi=find(Sr>srcrit(2) & Sr<=srcrit(3));
Rs(fi)=(kr(2)*mp(2))/(mum(2)-mp(2));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Zs(fi,1)=AA;
Bs(fi,2)=squeeze(Sr(fi))/mp(2);
Bs(fi,2)=Bs(fi,2)- (mz(2)*squeeze(Zs(fi,1)))/(ass(2)*mp(2)) - mz(2)/(ass(1)*grm(1));
Zs(fi,2)=0;
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));

% REGION 4
clear fi, fi=find(Sr>srcrit(3));
Bs(fi,1)=mz(1)/(ass(1)*grm(1));
Bs(fi,2)=mz(2)/(ass(2)*grm(2));
Zs(fi,1)=CC;
Zs(fi,2)=(ass(2)*Sr(fi))/mz(2)-mp(1)/grm(1)-mp(2)/grm(2) - CC;
Rs(fi)=(kr(2)*(grm(2)*Zs(fi,2)+mp(2)))./(mum(2)-(grm(2)*Zs(fi,2)+mp(2)));
Rscheck(fi)=(kr(1)*(grm(1)*Zs(fi,1)+mp(1)))./(mum(1)-(grm(1)*Zs(fi,1)+mp(1)));

%%%%%%%%%%%%%%%%%%
figure(1)
subplot(3,1,1)
plot(Sr,Rs,'b-','LineWidth',2);
hold on
plot(Sr,Rscheck,'r--','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end % for

subplot(3,1,2)
plot(Sr,squeeze(Bs(:,1)),'k-','LineWidth',2);
hold on
plot(Sr,squeeze(Bs(:,2)),'r-','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end 

subplot(3,1,3)
plot(Sr,squeeze(Zs(:,1)),'k-','LineWidth',2);
hold on
plot(Sr,squeeze(Zs(:,2)),'r-','LineWidth',2);
c=axis;
for ic=1:3
 plot([srcrit(ic) srcrit(ic)],[c(3) c(4)],'k:');
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% numerical soluiont %%%%%%%%%&
fisr=[3:6:150 ];
subplot(3,1,1), plot(Sr(fisr),Rs(fisr),'mo','LineWidth',2);
delt=[.1/24]; ntim=100000;
npp=0;
for is=[fisr];
  npp=npp+1;
  if (npp>4), npp=1; end
  srr=Sr(is);
  rsv(1)=1e-5; %Rs(is);
  for ip=1:nplank
   bsv(1,ip)=1e-5; %Bs(is,ip);
   zsv(1,ip)=1e-5; %Zs(is,ip);
  end % for ip
  tim(1)=0;
  for it=1:ntim
    rr=squeeze(rsv(it));
    bb=squeeze(bsv(it,:));
    zz=squeeze(zsv(it,:));
%
    gr(it)=0; gb(it,1:nplank)=0; gz(it,1:nplank)=0;
    for ip=1:nplank,
       gr(it)=gr(it)-mum(ip)*rr*bb(ip)/(rr+kr(ip));
       gb(it,ip)=gb(it,ip)+mum(ip)*rr*bb(ip)/(rr+kr(ip)) ...
                     -grm(ip)*zz(ip)*bb(ip)  ...
                     - mp(ip)*bb(ip);
       gz(it,ip)=gz(it,ip)+ass(ip)*grm(ip)*zz(ip)*bb(ip) ...
                     - mz(ip)*zz(ip);
    end % for ip 
    gr(it)=gr(it)+srr;
%%%
    rsv(it+1)=rr+gr(it)*delt;
    bsv(it+1,:)=bb+squeeze(gb(it,:))*delt; 
    zsv(it+1,:)=zz+squeeze(gz(it,:))*delt;
    tim(it+1)=tim(it)+delt;
  end % for it
  for ipl=1:2,
   figure(ipl)
   subplot(3,1,1), plot(Sr(is),rsv(it+1),'kx','LineWidth',2);
                   hold on
   subplot(3,1,2), plot(Sr(is),squeeze(mean(bsv(ntim-200:1:ntim,1))),'kx','LineWidth',2);
                   hold on
                  plot(Sr(is),squeeze(mean(bsv(ntim-200:1:ntim,2))),'rx','LineWidth',2);
   subplot(3,1,3), plot(Sr(is),squeeze(mean(zsv(ntim-200:1:ntim,1))),'kx','LineWidth',2);
                   hold on
                  plot(Sr(is),squeeze(mean(zsv(ntim-20:1:ntim,2))),'rx','LineWidth',2);

  end % for ipl
%%%%%%%
  if (npp==4),
   figure
   subplot(3,1,1)
   plot(tim, rsv,'k-','LineWidth',2);
   subplot(3,1,2)
   plot(tim, squeeze(bsv(:,1)),'k-','LineWidth',2);
   hold on
   plot(tim, squeeze(bsv(:,2)),'r--','LineWidth',2);
   subplot(3,1,3)
   plot(tim, squeeze(zsv(:,1)),'k-','LineWidth',2);
   hold on
   plot(tim, squeeze(zsv(:,2)),'r--','LineWidth',2);
  end % if npp
end % for is

