class NPZD extends Ecosystem {
  
  // NPZD model including light limitation; extension to DO is below.
  
  Statevar N,P,Z,D,par,totalN; 
  Flux Fuptake, Fingest, Fegest, Fpmort, Fzmort, Fexcret, Fremin;
  FreeParam mu0, ks, alpha, attSW, attP, m, I0, Ks, epsilon, fegest, zeta, Rremin;
  
  // system definition ----------------------
    
  NPZD() {
    N = addVar("nutrients","N","µM");
    P = addVar("phytoplankton","P","µM");
    Z = addVar("zooplankton","Z","µM");
    D = addVar("detritus","D","µM");
    totalN = addGhostVar("total nitrogen","T","µM");
    par = addGhostVar("PAR","PAR","W/m^2");
    
    Fuptake = addFlux("uptake","uptake", N, P);
    Fingest = addFlux("ingestion","ingest", P, Z);
    Fegest = addFlux("egestion","egest", P, D);
    Fexcret = addFlux("excretion","excret", Z, N);
    Fpmort = addFlux("P mortality","Pmort", P, D);
    Fzmort = addFlux("Z mortality","Zmort", Z, D);
    Fremin = addFlux("remineralization","Fremin", D, N);

    mu0 = addParam("P max growth","mu0","1/d",1.3,0,3);
    ks = addParam("nutrient half-saturation","ks","µM",1, 0.1, 5);
    alpha = addParam("initial P-E slope","alpha","", 0.15, 0.01, 0.2);
    attSW = addParam("light attenuation by seawater","attSW","1/m",0.1, 0, 0.2); // the Environment is repsonsible for using these to calculate the light profile
    attP = addParam("light attenuation by phytoplankton","attP","1/µM/m", 0.02, 0, 0.05);
    m = addParam("P mortality","m","1/d",0.1, 0, 0.4);
    I0 = addParam("Z max growth","I0","1/d", 4.8, 0, 6);
    Ks = addParam("grazing half-saturation","Ks","µM",3, 0.1, 10);
    Ks.logScale = true;
    epsilon = addParam("Z growth efficiency","epsilon","", 0.3, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.5, 0, 1);
    zeta = addParam("Z quadratic mortality","zeta","1/µM/d",2.0, 0, 4.0);
    Rremin = addParam("remineralization rate","remin","1/d",0.1, 0.01, 1);
    Rremin.logScale = true;
    // sinking gets set up by the Environment: don't turn it on here
  }
  
  
  void calcFluxes() {
    float small = 1e-6;
    for (int i=0; i<Ncells; i++) {
      float lightlim = alpha.current * par.current[i] / sqrt(small + sq(mu0.current) + sq(alpha.current * par.current[i]));
      float nutlim = N.current[i] / (small + ks.current + N.current[i]);
      Fuptake.current[i] = mu0.current * nutlim * lightlim * P.current[i];
      float Fgraz = I0.current * sq(P.current[i]) / (small + sq(Ks.current) + sq(P.current[i])) * Z.current[i];
      Fingest.current[i] = epsilon.current * Fgraz;
      Fegest.current[i] = (1 - epsilon.current) * fegest.current * Fgraz;
      Fexcret.current[i] = (1 - epsilon.current) * (1 - fegest.current) * Fgraz;
      Fingest.current[i] += Fexcret.current[i]; // because excretion is written as a flux from Z to N
      Fpmort.current[i] = m.current * P.current[i];
      Fzmort.current[i] = zeta.current * sq(Z.current[i]);
      Fremin.current[i] = Rremin.current * D.current[i];
      totalN.current[i] = N.current[i] + P.current[i] + Z.current[i] + D.current[i];      
    }
  }

}



// -----------------------------------------------------------------------------------



class NPZDO extends NPZD {

  Statevar Ox; 
  Flux Foxprod, Foxcons;
  FreeParam rOxNH4, rOxNO3, ks_ox_remin, ks_ox_zoop;
  
  // system definition ----------------------
  
  NPZDO() {
    Ox = addVar("oxygen","Ox","µM");
    Ox.min = 0;

    Foxprod = addFlux("oxygen production","oxprod", null, Ox);
    Foxcons = addFlux("oxygen consumption","oxcons", Ox, null);
    
    rOxNO3 = addParam("oxygen : nitrate ratio","rOxNO3","", 138./16., 0, 10);
    rOxNO3.hidden = true;
    rOxNH4 = addParam("oxygen : ammonium ratio","rOxNH4","", 106./16., 0, 10);
    rOxNH4.initial = rOxNO3.initial; // use rOxNO3 for everything, as if nitrification of NH4 is instantaneous
    rOxNH4.hidden = true;
    ks_ox_remin = addParam("oxygen half-saturation for remineralization","ks_ox_remin", "µM", 2, 0.1, 10); // small oxygen half-saturation for bacterial remineralization: Devol (1978) Deep Sea Res 25:137-146
    ks_ox_remin.hidden = true;
    ks_ox_zoop = addParam("oxygen half-saturation for zooplankton","ks_ox_remin", "µM", 30, 0.1, 100); // larger oxygen half-saturation for zooplankton metabolism, applied to all grazing fluxes
    ks_ox_zoop.hidden = true;
  }
    
  void calcFluxes() {
    super.calcFluxes();
    for (int i=0; i<Ncells; i++) {
      Foxprod.current[i] = rOxNO3.current * Fuptake.current[i];
      Foxcons.current[i] = rOxNH4.current * (Fexcret.current[i] + Fremin.current[i]);

      Fremin.current[i] *= Ox.current[i] / (Ox.current[i] + ks_ox_remin.current); 

      float grazSat = Ox.current[i] / (Ox.current[i] + ks_ox_zoop.current); 
      Fingest.current[i] *= grazSat;
      Fegest.current[i] *= grazSat;
      Fexcret.current[i] *= grazSat;
    }
  }

  
}



class NPZD_NetworkView extends NetworkView {
  
  NPZD eco;
  
  NPZD_NetworkView(Ecosystem eco) {
    this.eco = (NPZD)eco;
    varScale = 100;
    fluxScale = 4;
  }

  void define() {
    float theta0 = -120;
    addDot_polar(eco.D, 0.7, theta0-270, styles.medBlueColor);
    addDot_polar(eco.Z, 0.7, theta0-180, styles.dkBrownColor);
    addDot_polar(eco.P, 0.7, theta0-90, styles.medGreenColor);
    addDot_polar(eco.N, 0.7, theta0, styles.dkBlueColor);
    
    addArrow(eco.Fuptake, styles.ltGreenColor, 90);
    addArrow(eco.Fingest, styles.dkBrownColor, 90);
    addArrow(eco.Fegest, styles.medBlueColor, 30);
    addArrow(eco.Fexcret, styles.medBlueColor, 30);
    addArrow(eco.Fpmort, styles.medBlueColor, -30);
    addArrow(eco.Fzmort, styles.medBlueColor, 90);
    addArrow(eco.Fremin, styles.dkBlueColor, 90);
  }
  
}