class myNPZD extends dynamicalSystem {
  
  float totalN, totalN0; // internal variables
  int N_, P_, Z_, D_, par_, scale_, e0_, e1_;
  
  myNPZD() {
    define();
  }  
  
  dynamicalSystem clone() {
    return new myNPZD();
  }
  
  void define() {  
    adjustFluxesAtMaxMin = true;
    int numFramesToSave = 100;
    int maxvars = 11;
    int maxfluxes = 17;
    int maxparams = 40;
    allocate(numFramesToSave, maxvars, maxfluxes, maxparams);
    name = "my NPZD";
    shortname = "my NPZD";
    modelTimePerSecond = 2;
    timeUnits = "days";
    varUnits = "µM";
    tolerance = 1e-5;

    color Pcol = greenColor;
    color Zcol = redColor;
    color Dcol = color(60,70,130);
    color Ncol = color(10,15,70);
    
    float R = 0.8;
    N_ = addVar("nutrients","N", 0, Ncol, 0, -R);
    D_ = addVar("detritus","D",0, Dcol, R*sqrt(3)/2, -R);
    P_ = addVar("phytoplankton","P", 0, Pcol, -R*sqrt(3)/2, 0.5*R);
    Z_ = addVar("zooplankton","Z", 0, Zcol, R*sqrt(3)/2, 0.5*R);
    par_ = addVar("PAR","PAR", 100, yellowColor, 0,0); // this is where communicationAfterCalcStep() stores PAR for each cell
    vars[par_].display = false;
    e0_ = addEndpoint("N source pool", "", 0, -1.3*R);
    e1_ = addEndpoint("D oblivion", "",  R*sqrt(3)/2, -1.3*R);
    
    // scale bar: off by default
    color scaleColor = transparency(color(255),0.5);
    scale_ = addVar("scale","1", 1, color(255), -2*R,R*2);
    vars[scale_].display = false;
    vars[scale_].constant = true;
    scale_ = addVar("scale","5", 5, color(255), -0.6*R,R*2);
    vars[scale_].display = false;
    vars[scale_].constant = true;
    scale_ = addVar("scale","10", 10, color(255), R,R*2);
    vars[scale_].display = false;
    vars[scale_].constant = true;
    // note: the last var associated with scale_ is the one that actually sets the scale for drawing later on
    
    setInitialConditions();
        
    addFlux("N supply","Nsupply", e0_, N_, Ncol, 0);     
    addFlux("N uptake","upt", N_, P_, Pcol,120);
    addFlux("P mortality","Pmort", P_, D_, Dcol, 60);
    addFlux("ingestion","ingest", P_, Z_, Zcol,120);
    addFlux("egestion","egest", P_, D_, Zcol,120);
    addFlux("basal metabolism","metab", Z_, N_, Ncol,110);
    addFlux("ingestion-related excretion","excret", Z_, N_, Zcol,120);
    addFlux("Z mortality","Zmort", Z_, D_, Dcol,60);
    addFlux("remineralization","remin", D_, N_, Ncol,20);
    addFlux("D removal","Dremoval", D_, e1_, Dcol, 0);

    addParam("attenuation by seawater (1/m)","attSW", 0.13, 0, 0.13, Pcol);
    addParam("attenuation by chl (1/m/(mg/m3))","attChl", 0.0072, 0, 0.02, Pcol);
    addConstant("PAR fraction","parFrac", 0.43);
    addConstant("chl:N ratio","chl2N", 2.5); // this is handled dynamically in roms
    addConstant("growth-light initial slope (alpha)","alpha",0.07); // in roms, PhyIS
    addParam("N injection/removal rate (µM/day)", "S", 0, -2, 2, Ncol);
    addParam("P max 24 h growth rate (1/day)","Vm",1.3, 0,3, Pcol);
    addParam("half-saturation for N uptake (µM)","k",4.6, 0.1, 5, Pcol); // the inverse of this appears in bioFasham.h
    addParam("P mortality rate (1/day)","mp",0.1, 0, 0.2, Pcol);
    addParam("Z max grazing rate (1/day)","Imax",4.8, 0,10, Zcol);
    addParam("half-saturation for grazing (µM)","K",3, 0.1, 10, Zcol); // in roms, the SQUARE of this is specified (...P^2/(K+P^2)...)
    addParam("egestion fraction","gamma",0.35, 0, 1, Zcol); // in roms, 1-ZooAE_N
    addParam("excretion fraction ","excr",0.35, 0, 1, Zcol);
    addConstant("Z basal metabolism (1/day)","basal",0);
    addParam("Z mortality rate (1/µM/day)","mz",2, 0, 6, Zcol);
    addParam("D remineralization rate (1/day)","rem",0.1, 0, 1, Dcol);
    addParam("D sinking (m/day)","wsinkD",8, 0, 20, Dcol);
    
    RESPOND_TO_PARAMCHANGE();
    
  }
  
  void setInitialConditions() {
    setInitialConditions(10);
  }
  
  void setInitialConditions(float romsNO3) {
    vars[P_].initial = 0.1;
    vars[Z_].initial = 0.01;
    vars[N_].initial = romsNO3;
    vars[D_].initial = 0;
  }
  
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    
    float P = vslice[P_];
    float Z = vslice[Z_];
    float N = vslice[N_];
    float D = vslice[D_];
    float PAR = vslice[par_];
    
    // N injection & flowthrough
    F[find("Nsupply")] = readout("S");
    if (systems[currSysID] == this) { // if this is _not_ a child of another system
      F[find("Dremoval")] = readout("wsinkD") / 20;
    } else {
      F[find("Dremoval")] = 0;
    }
    
    // phyto growth
    float alphaI = readout("alpha") * PAR + tolerance;
    float Vm = readout("Vm");
    float mu = Vm * alphaI / sqrt(sq(Vm) + sq(alphaI));
    float k = readout("k") + tolerance;
    F[find("upt")] = mu * N / (k + N) * P;
    
    // grazing
    String grazingResponse = "quadratic";
    float K = readout("K") + tolerance; // grazing half-sat for Z and Z2
    float grazing = 0;
    float grazing2 = 0;
    if (grazingResponse.equals("linear")) {
      grazing = readout("Imax") * P / (K + P) * Z;
    } else if (grazingResponse.equals("quadratic")) {
      grazing = readout("Imax") * sq(P) / (sq(K) + sq(P)) * Z;
    }
    float gamma = readout("gamma");
    float excrFrac = readout("excr");
    F[find("ingest")] = (1-gamma) * grazing;
    F[find("egest")] = gamma * grazing;
    F[find("excret")] = excrFrac * grazing;
    
    // other loss terms
    F[find("metab")] = readout("basal") * Z;
    float mp = readout("mp");
    F[find("Pmort")] = mp * P;
    float mz = readout("mz");
    F[find("Zmort")] = mz * sq(Z);
    F[find("remin")] = readout("rem") * D;
    
    return F; 
  }
  
  void adjustments_after_calc_step() {
    totalN = vars[N_].current + vars[P_].current + vars[Z_].current + vars[D_].current;
    for (int i=0; i<nvars; i++) {
      if (i != par_) {
        vars[i].min = 0.001;
// the *100 is a kluge to solve a now-forgotten problem. Maybe even a now-solved problem.
        vars[i].max = totalN*100;
        vars[i].scale = 5;
      }
    }
  }
  
  void RESPOND_TO_PARAMCHANGE() {
    vars[D_].sinkingRate = readout("wsinkD");
  }
  
  float totalP() {
    return vars[P_].current;
  }
  
  float totalNitrogen() {
    return totalN;
  }
  
  void showDiagnostics() {
    float Fup = fluxes[find("upt")].current;
    float mu = Fup / totalP();
    float Fgraz = fluxes[find("ingest")].current + fluxes[find("egest")].current;
    float g = Fgraz / totalP();
    setFont("flexi");
    textSize(36);
    textLeading(40);
    textAlign(CENTER);
//    String diagnosticsString = "P = " + str_sigfigs(totalP()) + " µmol/L" + "\n" + "µ = " + str_sigfigs(mu) + "/day" + "\n" + "g = " + str_sigfigs(g) + "/day";
    String diagnosticsString = "chl Å " + str_sigfigs(2.93*totalP()) + "\n" + "µ = " + str_sigfigs(mu) + "\n" + "g = " + str_sigfigs(g);
    // N:chl ratio of 2.93 from Raphe nov 2006, for RISE3
    
    // blank out a space for the text
    fill(transparency(color(255),0.6));
    noStroke();
    float dx = 1.1 * textWidth(diagnosticsString);
    float dy = 3 * 40 * 1.1;
    rectMode(CENTER);
    rect(networkRect[0]+networkRect[2]/2, networkRect[1]+networkRect[3]/2, dx, dy);
    rectMode(CORNER);
    
    // redraw the var dots to highlight them
    draw_network();
    
    // write the text
    setFont("flexi");
    textSize(36);
    textLeading(40);
    textAlign(CENTER);
    fill(transparency(color(255),0.8));
    text(diagnosticsString,networkRect[0]+networkRect[2]/2,networkRect[1]+networkRect[3]/2-28);
    fill(transparency(orangeColor,0.8));
    text(diagnosticsString,networkRect[0]+networkRect[2]/2,networkRect[1]+networkRect[3]/2-28);
  }
  
  void adjustments_before_drawing() { // set dot positions so they cluster together
    if (systems[currSysID] != this) { // if this is a child of another system (as opposed to a 0-D model run on its own)
      vars[P_].xc = -vars[P_].radius;
      vars[P_].yc =  vars[P_].radius;
      vars[Z_].xc =  vars[Z_].radius;
      vars[Z_].yc =  vars[Z_].radius;
      vars[N_].xc = -vars[N_].radius;
      vars[N_].yc = -vars[N_].radius;
      vars[D_].xc =  vars[D_].radius;
      vars[D_].yc = -vars[D_].radius;    
    }
  }

} // myNPZD