clear all
close all

toppmin=-10;
toppmax=45;

tlocal=toppmin:.1:toppmax;
trange=toppmax-toppmin;

% pseudo-Arrenhius
     NEW_A=1; % need to chose
     NEW_AE=-4000;
     NEW_Tref=293.15;
     NEW_K=273.15;
     NEW_C= 4.d0;
     NEW_tau2 = 0.0;
     NEW_tau1 = 0.800 % 1.d0/(1*1.7);

figure, subplot(2,2,1), hold off
for isize=1:1
 if (isize==1), NEW_B=3e-4; cstr=['r-']; end
 if (isize==2), NEW_B=1e-4; cstr=['b--']; end
 for ipt=1:10,
   phytoTempOptimum=-3+(ipt-1)*4.;
   sphytoTempOptimum(ipt)=phytoTempOptimum;
%  phytoTempOptimum=7+ipt*.5;

     for it=1:length(tlocal),
            NEW_arren(it)=(NEW_A*exp(NEW_AE*(1/(tlocal(it)+NEW_K)-1/NEW_Tref)));
            NEW_tlimit(it) =( NEW_arren(it)* ...
             exp(-NEW_B* ...
          (abs(tlocal(it) - phytoTempOptimum))^NEW_C));
            NEW_tlimit(it) = NEW_tlimit(it) -NEW_tau2;
            NEW_tlimit(it) = NEW_tau1*max(NEW_tlimit(it), 1d-10);
            NORM_arren(it) = NEW_arren(it)  -NEW_tau2;
            NORM_arren(it) = NEW_tau1*max(NORM_arren(it), 1d-10);
     end % for it
%  plot(tlocal,NEW_tlimit,cstr,'linewidth',2.); hold on
   plot(tlocal,NORM_arren,'k-','linewidth',2.); hold on
   plot([-2 40],[1 1],'k:');
   axis([-2 40 0 1.4]);
   bolder
   ylabel('growth factor (unitless)');
   xlabel('temperature (^oC)');
 end % for ipt
end % for isize

no3=0.005:0.005:0.5;
kn=8.2645515845713290e-3;
mumax=1;
nlimit=no3./(no3+kn);
subplot(2,2,2),
plot(no3,nlimit,'k-','linewidth',2.);
axis([0 0.5 0 1.4])
   bolder
   ylabel('growth factor (unitless)');
   xlabel('NO3 (mmol/m3)');


for it=1:length(tlocal),
 for inn=1:length(no3);
   totlimit(it,inn)=NORM_arren(it)*nlimit(inn);
 end %
end %


fid=fopen('/home/stephdut/Input/input_ecco/lev01_nitrate_ann-3d.bin','r','b');
nit=fread(fid,'float32');
nit=reshape(nit,360,160,23);
   istart=78936;
   timestep=24*30; timestep1=24;

      nx=360; ny=160; nz=23;
   theta=zeros(nx,ny,nz); 
   for itt=1:360,
    time1=istart+(itt-1)*timestep1;
    if (time1<1000000); timestr=['0000',num2str(time1)]; end
    if (time1<100000); timestr=['00000',num2str(time1)]; end
    if (time1<10000); timestr=['000000',num2str(time1)]; end
    if (time1<1000); timestr=['0000000',num2str(time1)]; end
    if (time1<100); timestr=['00000000',num2str(time1)]; end
%
    fid=fopen(['/nobackup1/stephdut/ECCO_Iter3.73daily/DDtheta/DDtheta.',timestr,'.data'],'r','b');
    field=fread(fid,'float32'); fclose(fid);
    field=reshape(field,nx,ny,nz);
    theta=theta+field/360;
   end % for itt


subplot(2,2,3)

x1=180; x2=220; y1=105; y2=130;
stlocal=tlocal(250:345);
snlocal=fliplr(no3(5:length(no3)));
clear tmp, tmp=squeeze(theta( x1:x2, y1:y2 ,1));
clear fi, fi=find(tmp>0); stheta=tmp(fi);
clear tmp,tmp=squeeze(nit( x1:x2, y1:y2 ,1)); snit=tmp(fi);

plot(stheta,snit,'k.')
hold on
plot(stlocal,snlocal,'linewidth',2.);
bolder
axis([15 25 0 0.5]);
   xlabel('temperature (oC)');
   ylabel('NO3(mmol/m3)');


%%%%%%


for it=1:length(stlocal),
   stotlimit(it)=NORM_arren(249+it)*(snlocal(it)/(snlocal(it)+kn));
end

subplot(2,2,4),
plot(stlocal,stotlimit,'k-','linewidth',2.); 
bolder
axis([15 25 0 1]);
   xlabel('temperature (oC)');
   ylabel('combined growth factor');

%
figure
colormap default, cmap=colormap; colormap(cmap([18:2:57],:));
pcolor(tlocal,no3,totlimit'); shading flat, colorbar
axis([-2 40 0 0.5]);
   bolder
   xlabel('temperature (oC)');
   ylabel('NO3 (mmol/m3)');
   title('combined growth factor');


figure
colormap default, cmap=colormap; colormap(cmap([18:2:57],:));
subplot(2,2,1)
pcolor(squeeze(theta(x1:x2,y1:y2,1))'); shading flat, colorbar
subplot(2,2,2)
pcolor(squeeze(nit(x1:x2,y1:y2,1))'); shading flat, colorbar

