#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.014 + (0.044-0.014) * 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           absorbence at 490nm
      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 ration */
  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;
}