#include <math.h>

/*--------------------------------------------------------------------------*/

double opp_chlCarbon2( double aph,
		       double bbp,
		       double irr,
		       double k490,
		       double mld,
		       double daylength) {
/*

   !Description:     opp_chlCarbon - computes daily primary productivity using
                     a chl:Carbon ratio.  The ChlCarbon algorithm estimates productivity 
		     using aph (m-1), bbp (m-1), surface irradiance (Einsteins m-2 d-1),
                     k490 (m-1), mld (m), and day length (hours).

growth rate (u) = u(max) * f(nut,T) * f(Ig)

where:

u(max) = 2

f(nut,T) = (Chl/C)sat / (Chl/C)max

     where: (Chl/C)sat = SeaWiFS Chl:C = converted ratio of aph/bbp
               (Chl/C)max = 0.022 + (0.045-0.022) * exp (-3 * Ig)

         where: Ig = median mixed layer light level = surface irradiance * 
	                                              exp (-k490 * MLD/2)

f(Ig) = 1 - exp (-3 * Ig)



   !Input Parameters:  
      aph            absorption due to phytoplankton
      bbp            backscatter
      irr            Photosynthetically available radiation in Einsteins per
                     day per square meter
      k490           diffuse attenuation coefficient at 490 nm (units of 1 / m )
      mld            mixing layer depth in meters
      dayL           Length day in decimal hours.

   !Output Parameters: 
             Primary productivity in milligrams Carbon per square meter
                     per day

   !Revision History:  

      June-16-2005 first version (Robert O'Malley)
      
      After program opp_chlCarbon.c
                and opp_befa.c
		 
        replace the z_eu calc in opp_chlCarbon.c
	with the method used in opp_befa.c

   !References and Credits
   
      Behrenfeld,M.J; Boss,E.; Siegel,D.A.; Shea,D.M.; 2005. Carbon-based Ocean
      Productivity and Phytoplankton Physiology from Space.  Global Biogeochemical
      Cycles, Volume 19
      
*/

  double uMax;			/* max growth rate */
  double u;			/* growth rate */
  double chlCarbonMax;		/* max chl:carbon ratio */
  double nutTempFunc;		/* f(nut,T) */
  double chlCarbonSat;		/* satellite chl:carbon ratio */
  double chl;			/* satellite chl */
  double carbon;		/* bbp converted to carbon */
  double Ig;			/* median mixed layer light level */
  double IgFunc;		/* f(Ig) */
  double irrFunc;		/* irradiance function describing water column */
  double z_eu;			/* euphotic depth at 1% light level */
  double npp;
  double chl_tot;
  
  if(irr <= 0.0)
    return 0.0;

  uMax = 2.0;
  
  Ig = irr / daylength * exp(-k490 * mld / 2.0);

  IgFunc = 1 - exp(-3.0 * Ig);

  chl = aph / 0.05582;

  if(bbp < 0.00035)
    bbp = 0.00036;
  carbon = 13000.0 * (bbp - 0.00035);
  
  chlCarbonSat = chl / carbon;

  chlCarbonMax = 0.022 + (0.045-0.022) * exp(-3.0 * Ig);
  
  nutTempFunc = chlCarbonSat / chlCarbonMax;

  /* irr function */
  irrFunc = 0.66125L * irr / ( irr + 4.1L );
  
  u = uMax * nutTempFunc * IgFunc;
  
   /* Calculate euphotic depth (z_eu) with Morel's Case I model.            */
   /* Calculate chl_tot from Satellite Surface Chlorophyll Data.            */

  if (chl <  1.0L) 
    chl_tot = 38.0L * pow( chl, 0.425 );
  else
    chl_tot = 40.2L * pow( chl, 0.507 );

   
  z_eu = 200.0L * pow( chl_tot, (-.293) );

  if (z_eu <= 102.0L) 
    z_eu = 568.2L * pow( chl_tot, (-.746) );

  //  z_eu = -log(0.01) / k490;

  npp = carbon * u * irrFunc * z_eu;

  return npp;
}