%set of differential equations of
%oschlies & garcon 1999  npzd model

global w no3e sra  sst lat K Jave 

sda = 60*60*24; %seconds in a day
sda = 1 ;  % per day
sst = 10.;
lat=-50;
%K=1;
MLD=25;  %1;  %25;
%sra = 300;

%************* parameter list *************
%pyhtoplankton coefficients
 alpha = 0.025/sda;            %Phyto: [1/s/(W/m^2)] Initial slope of P-I curve
 par   = 0.43;                 %Phyto: [] Photosynthetically active radiation
 kw    = 25;                   %Phyto: [m] Light attenuation due to water
 kc    = 0.03;                 %Phyto: [1/m(mmolN/m^3)^-1] Light att. due to Phy.
 a     = 0.6/sda;              %Phyto: [1/s] Max. growth rate parameter
 b     = 1.066;                %Phyto: [] Max. growth rate parameter
 c     = 1;                    %Phyto: [1/(degC)] Max. growth rate parameter
 K1    = 0.5;                  %Phyto: [mmol/m^3] Half-sat. const. for N uptake
 muep  = 0.03/sda;             %Phyto: [1/s] Specific mortality rate
% muepst= 0.025/sda;           %Phyto: [(mmolNm^-3)^-1/s] Quad. mort. rate
 Jmax  = a*b^(c*sst);          %Phyto: [1/s] light-saturated growth oa+gar 1999
 %zooplankton coefficients
  gam1  = 0.75;                %Zoo: [] Assim. efficiency
  g     = 2/sda;               %Zoo: [1/s] Max. graz. rate
  eps   = 1/sda;               %Zoo: [(mmolNm^-3)^-2/s] Prey capture rate
  muez  = 0.2/sda;             %Zoo: [mmolm^-3)^-1/s] Quad. mortality
  gam2  = 0.03/sda;            %Zoo: [1/s] Excretion
%detrial coefficients
  mued  = 0.05/sda;            %Det: [1/s] Remineralization rate
  w_sink= 5/sda;              %Det: [m/s] Sinking velocity 
%*****************************************


%////////calculation of a couple of terms///////////
%light calculation
       fx1    = 2.0*pi*K./12;
       declin = 0.006918     ...
             -0.399912*cos(   fx1)+0.070257*sin(   fx1)    ... 
             -0.006758*cos(2.*fx1)+0.000907*sin(2.*fx1)    ...    
             -0.002697*cos(3.*fx1)+0.001480*sin(3.*fx1);   
       fx1    = pi/180. * lat; 
       cobeta = max([0. sin(fx1)*sin(declin)+cos(fx1)*cos(declin)]);
       cobeta = sqrt(1.-(1.-cobeta.^2.)/1.33.^2.);
       fx2    = max([-1.  min( 1.,-tan(fx1)*tan(declin))] );
       daylen = max([1.e-12 acos(fx2)/pi]);
       rayb(1)= 1.0;

       rayb(2)= max([1e-15 exp(-MLD/(cobeta*kw))]);

          fx1 = 0.5*daylen*daylen/alpha;  
	   sra = 2*sra/daylen;
          radbio = max([1.0 par*sra]);
          vpbio = a*b.^(c*sst);
          fx3 = fx1*vpbio/(radbio*rayb(1));
          fx4 = fx1*vpbio/(radbio*rayb(2));
          fu1 = sqrt(fx3.^2.+daylen.^2.);
          fu1 = fu1-daylen*log((daylen+fu1)/fx3);
          fu2 = sqrt(fx4.^2.+daylen.^2.);
          fu2 = fu2-daylen*log((daylen+fu2)/fx4);
       Jave = -1.*vpbio/MLD * cobeta*kw  ... 
                 *(fu1 - fu2 - fx3 + fx4);
	Jave

