#include <math.h>
/*--------------------------------------------------------------------------*/
double opp_cbpm2( double chl,
double bbp,
double irr,
double k490,
double mld,
double zno3,
double daylength) {
/*
!Description: opp_cbpm2 - computes daily primary productivity using a chl:Carbon ratio.
This is a spectrally resolved version of the cbpm, using nine separate
wavelengths. It is also depth resolved, integrating the effects from
the surface down to a fixed depth of 300 m.
The cbpm2 algorithm estimates productivity using chl (m-1), bbp (m-1),
surface irradiance (Einsteins m-2 d-1), k490 (m-1), mld (m), zno3 (m)
and day length (hours).
Net primary productivity is carbon * growth rate, where carbon is proportional to particulate
backscatter
carbon = 13000 * (bbp - 0.00035)
and growth rate is a function of nutrient and temperature stress (f(nut,T) and photoacclimation
(f(Ig))
growth rate (u) = umax * f(nut,T) * f(Ig)
where:
umax = 2
f(nut,T) = ((Chl/C)sat - y0) / ((Chl/C)max - y0)
f(Ig) = 1 - exp (-5 * Ig)
and:
(Chl/C)sat = ratio of satellite observed chl and carbon (carbon from bbp)
(Chl/C)max = 0.022 + (0.045-0.022) * exp (-3 * Ig)
Ig = median mixed layer light level
= surface irradiance * exp (-k(lambda) * MLD/2)
The above items are analyzed for nine separate wavelengths, and is vertically resolved to a depth
of 300 m.
For more details, please see the paper by Westberry, et al (2008)
!Input Parameters:
chl chlorophyll concentration
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
zno3 depth of the nitrocline
daylength length of the day in decimal hours.
!Output Parameters:
Primary productivity in milligrams Carbon per square meter
per day
!Dependencies:
function austinPetzold_1986 ( double lambda, double K490 )
given a reference k490 vlaue, determine k(lambda) for a specified lambda
ref:
Austin, R. W., and T. J. Petzold (1986), Spectral dependence of the diffuse
attenuation coefficient of light in ocean waters, Opt. Eng., 25, 473 – 479
!Revision History:
08-16-2010 first release version (Robert O'Malley)
[original code written in matlab by T. Westberry]
01-05-2011 O'Malley
add uMax trap on mu[m]
correct z_eu determination
01-25-2023 O'Malley
adjusted the depth of integration
from 200 to 300 m
per request by T. Westberry
!References and Credits
Westberry, T. Behrenfeld, M.J., Siegel, D.A., and Boss, E.; 2008. Carbon-based
primary productivity modeling with vertically resolved photoacclimation. Global
Biogeochemical Cycles, Vol. 22, GB2024, doi:10.1029/2007GB003078
*/
double austinPetzold_1986( double, double );
double uMax; /* max growth rate */
double chlCarbonMax; /* max chl:carbon ration */
double nutTempFunc; /* f(nut,T) */
double chlCarbonSat; /* satalite chl:carbon ratio */
double carbon; /* bbp converted to carbon */
double IgFunc; /* f(Ig) */
double IgFuncz; /* f(Ig) below the mixed layer depth */
double z_eu; /* euphotic depth at 1% light level */
double npp; /* net primary production */
/* --------------------- */
/* spectral variables */
/* --------------------- */
double lambda[] = { 400, 412, 443, 490, 510, 555, 625, 670, 700 };
double parFraction[] = { 0.0029, 0.0032, 0.0035, 0.0037, 0.0037, 0.0036, 0.0032, 0.0030, 0.0024 };
double X[] = { .11748, .122858, .107212, .07242, .05943, .03996, .04000, .05150, .03000 };
double e[] = { .64358, .653270, .673358, .68955, .68567, .64204, .64700, .69500, .60000 };
double Kw[]= { .01042, .007932, .009480, .01660, .03385, .06053, .28400, .43946, .62438 };
double Kd[9];
double Kbio;
double Kdif[9];
double Klambda[9];
double Eo[9];
double Ez_mld[9];
double par_mld;
double delChlC;
double y0;
/* --------------------------- */
/* depth resolved variables */
/* --------------------------- */
double z[300]; /* depths */
double chl_C[300]; /* chl:c ratio */
double chlz[300]; /* chl */
double mu[300]; /* growth */
double Ezlambda[9][300]; /* fraction of light at nine wavelengths */
double parz[300]; /* total light */
double prcnt[300]; /* percent light */
double Cz[300]; /* carbon */
double ppz[300]; /* npp */
int i;
int m;
int mzeu;
double r;
double prcnt0;
double prcnt1;
double z0;
double z1;
double numerator;
double denominator;
double fraction;
double deltaZ;
if(irr <= 0.0){
return 0.0;
}
/* --------------------- */
/* initialize values */
/* --------------------- */
z_eu = -9999; // 1.05.2011
y0 = 0.0003; /* min chl:c when mu = 0 */
for (i = 0; i<300; i++){
z[i]= (float)(i+1);
}
r = 0.1;
uMax = 2.0; /* after BANSE (1991) */
npp = 0.0;
mzeu = 0;
for (i=0; i<9; i++) {
Klambda[i]= austinPetzold_1986(lambda[i],k490);
Eo[i] = irr * parFraction[i];
Ez_mld[i] = Eo[i] * 0.975 * exp(-Klambda[i] * mld / 2.0);
}
/* ----------------------------- */
/* reintegrate to get par at */
/* depth ... */
/* do trapezoidal integration */
/* ----------------------------- */
par_mld = 0.0;
for (i = 0; i < 8; i++ ){
par_mld += (lambda[i+1]-lambda[i])*(Ez_mld[i+1]+Ez_mld[i])/2;
}
par_mld /= daylength;
IgFunc = 1 - exp(-5.0 * par_mld);
if(bbp < 0.00035)
bbp = 0.00036;
carbon = 13000.0 * (bbp - 0.00035);
chlCarbonSat = chl / carbon;
if ( chlCarbonSat < y0 ) {
chlCarbonSat = y0;
}
chlCarbonMax = 0.022 + (0.045-0.022) * exp(-3.0 * par_mld);
delChlC = chlCarbonMax - chlCarbonSat;
nutTempFunc = (chlCarbonSat - y0) / (chlCarbonMax - y0);
/* ''''''''''''''''''''''''' */
/* calculate Kd offset */
/* carry through to depth */
/* non-chl attenuation */
/* ------------------------- */
for (i=0; i<9; i++) {
Kbio = X[i] * pow(chl,e[i]);
Kd[i] = Kw[i] + Kbio;
Kdif[i] = Klambda[i] - Kd[i];
}
/* ''''''''''''''''''''''''''''''''''' */
/* integrate down the water column */
/* in one-meter steps */
/* ----------------------------------- */
for (m=0; m<300; m++) {
/* ---------------------------------------------- */
/* if you are in the mixed layer, do this way */
/* ---------------------------------------------- */
if ( z[m] < mld ) {
chl_C[m] = chlCarbonSat;
chlz[m] = chl_C[m] * carbon;
mu[m] = uMax * nutTempFunc * IgFunc;
if ( mu[m] > uMax ) { // 1.05.2011
mu[m] = uMax; // 1.05.2011
} // 1.05.2011
for ( i=0; i<9; i++) {
Ezlambda[i][m] = Eo[i]*0.975*exp(-Klambda[i]*z[m]);
}
parz[m] = 0.0;
for (i = 0; i < 8; i++ ){
parz[m] += (lambda[i+1]-lambda[i])*(Ezlambda[i+1][m]+Ezlambda[i][m])/2;
}
Cz[m] = carbon;
} else {
/* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''' */
/* if below mixed layer must treat properties differently */
/* ---------------------------------------------------------- */
for (i=0; i<9; i++) {
Kbio = X[i] * pow(chlz[m-1],e[i]); /* after Morel & Maritorena (2001) */
Kd[i] = Kw[i] + Kbio;
Kd[i] += Kdif[i];
Ezlambda[i][m] = Ezlambda[i][m-1]*exp(-Kd[i]*1.0);
}
parz[m] = 0.0;
for (i = 0; i < 8; i++ ){
parz[m] += (lambda[i+1]-lambda[i])*(Ezlambda[i+1][m]+Ezlambda[i][m])/2;
}
deltaZ = zno3 - z[m];
if ( deltaZ < 0 ) {
deltaZ = 0;
}
chl_C[m] = (0.022 + (0.045-0.022) * exp(-3.0 * parz[m] / daylength));
chl_C[m] -= delChlC * (1-exp(-0.075*deltaZ));
IgFuncz = 1 - exp(-5.0 * parz[m]/daylength);
mu[m] = uMax * nutTempFunc * IgFuncz;
if ( mu[m] > uMax ) { // 1.05.2011
mu[m] = uMax; // 1.05.2011
} // 1.05.2011
if (mu[m-1] >= r ) {
Cz[m] = carbon;
} else {
Cz[m] = carbon * mu[m-1] / r;
}
chlz[m] = chl_C[m] * Cz[m];
}
prcnt[m] = parz[m] / (irr * 0.975);
/* track this to get to the euphotic depth */
if ( prcnt[m] >= 0.01 ) {
mzeu = m;
} else {
/* ''''''''''''''''''''''''''' */
/* now find 1% light depth */
/* in case the user wants */
/* to use this information */
/* --------------------------- */
if (z_eu == -9999 ) { // 01.05.11
prcnt0 = prcnt[mzeu];
prcnt1 = prcnt[mzeu+1];
z0 = z[mzeu];
z1 = z[mzeu+1];
numerator = prcnt0 - 0.01;
denominator = prcnt0 - prcnt1;
fraction = numerator / denominator;
z_eu = z0 + (z1-z0)*fraction;
}
}
ppz[m] = mu[m] * Cz[m];
}
/* ------------------------------- */
/* do trapezoidal integration */
/* from m = 0 to m = 300 */
/* ------------------------------- */
// note: 186 m is the euphotic depth for pure water
if ( mzeu < 186 ) {
npp = 0;
for (i = 0; i < 299; i++ ){
npp += (z[i+1]-z[i])*(ppz[i+1]+ppz[i])/2;
}
} else {
npp = -9999;
}
return npp;
}
/* ================================================================= */
double austinPetzold_1986 ( double lambda,
double K490 ) {
double wave[] = { 350, 360, 370, 380, 390, 400,
410, 420, 430, 440, 450, 460, 470, 480, 490, 500,
510, 520, 530, 540, 550, 560, 570, 580, 590, 600,
610, 620, 630, 640, 650, 660, 670, 680, 690, 700 };
double M[] = { 2.1442, 2.0504, 1.9610, 1.8772, 1.8009, 1.7383,
1.7591, 1.6974, 1.6108, 1.5169, 1.4158, 1.3077, 1.1982, 1.0955, 1.0000, 0.9118,
0.8310, 0.7578, 0.6924, 0.6350, 0.5860, 0.5457, 0.5146, 0.4935, 0.4840, 0.4903,
0.5090, 0.5380, 0.6231, 0.7001, 0.7300, 0.7301, 0.7008, 0.6245, 0.4901, 0.2891 };
double Kdw[] = { 0.0510, 0.0405, 0.0331, 0.0278, 0.0242, 0.0217,
0.0200, 0.0189, 0.0182, 0.0178, 0.0176, 0.0176, 0.0179, 0.0193, 0.0224, 0.0280,
0.0369, 0.0498, 0.0526, 0.0577, 0.0640, 0.0723, 0.0842, 0.1065, 0.1578, 0.2409,
0.2892, 0.3124, 0.3296, 0.3290, 0.3559, 0.4105, 0.4278, 0.4521, 0.5116, 0.6514 };
double l0;
double l1;
double k0;
double k1;
double m0;
double m1;
double kdiff;
double mdiff;
double num;
double den;
double frac;
double Kdw_l;
double M_l;
double Kd;
int ref;
int i;
// -- INTERPOLATE TO WAVELENGTH OF INTEREST -- //
for (i = 1; i < 36; i++) {
if ( wave[i] >= lambda ) {
l1 = wave[i];
k1 = Kdw[i];
m1 = M[i];
l0 = wave[i-1];
k0 = Kdw[i-1];
m0 = M[i-1];
break;
}
}
num = lambda - l0;
den = l1 - l0;
frac = num / den;
kdiff = k1 - k0;
Kdw_l = k0 + frac*kdiff;
mdiff = m1 - m0;
M_l = m0 + frac*mdiff;
// -- GET REFERENCE WAVELENGTH (=490 FOR NOW) AND APPLY MODEL -- //
ref = 14;
Kd = (M_l/M[ref]) * (K490 - Kdw[ref]) + Kdw_l;
return Kd;
}
|
