class NPZsystem extends dynamicalSystem {
  
  float totalN, totalN0; // internal variables
  int N_, P_, Z_;
  float ln2 = 0.6931;
  
  NPZsystem() { // separating the constructor and define() means you can easily make a variant that calls the same define() but does other things in the constructor
    define();
  }  
    
  void define() {  
    int numFramesToSave = 100;
    int maxvars = 3;
    int maxfluxes = 12;
    int maxparams = 14;
    allocate(numFramesToSave, maxvars, maxfluxes, maxparams);
    name = "NPZ";
    shortname = "NPZ";
    totalN0 = 8;
    float P0 = 1; float Z0 = 1;
    modelTimePerSecond = 4;
    timeUnits = "days";
    varUnits = "µmol nitrogen";
    tolerance = 1e-5;

    color Pcol = greenColor; color Zcol = redColor;
    color Ncol = color(20,30,90);
    float R = 1; float po = PI/180;
    N_ = addVar("nutrients","N", totalN0-P0-Z0, Ncol, 0,-R);
    P_ = addVar("phytoplankton","P", P0, Pcol,R*cos(150*po),R*sin(150*po));
    Z_ = addVar("zooplankton","Z", Z0, Zcol,R*cos(30*po),R*sin(30*po));
    
    addFlux("uptake","upt", N_, P_, Pcol,120-20);
    addFlux("ingestion","ingest", P_, Z_, Zcol,120);
    addFlux("P mortality","Pmort", P_, N_, Ncol,120);
    addFlux("Z baseline losses","Zmort", Z_, N_, Ncol,120);
    addFlux("Z ingestion-related losses","excret", Z_, N_, darken(Pcol),120-30);
    
    addParam("P max growth (doublings/day)","Vm",1,0,1.5, Pcol);
    addParam("Z max growth (doublings/day)","Rgr",0.75, 0, 1.5, Zcol);
    addParam("nutrient uptake half-saturation (µmol/L)","k",5, 0.1, 5, Pcol);
    addParam("P mortality (1/day)","m",0.1, 0, 0.4, Pcol);
    addParam("grazing half-saturation (µmol/L)","K",2, 0, 4, Zcol);
    addParam("fraction of prey re-excreted","Rex",0.5, 0, 1, Zcol);
    addParam("Z baseline losses (1/day)","g",0.2, 0, 0.4, Zcol);
  }
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    float K = readout("K"); float k = readout("k"); float Vm = readout("Vm");
    float m = readout("m");
    
    F[find("upt")] = Vm * ln2 * vslice[N_]/(k+vslice[N_]) * vslice[P_];
    F[find("ingest")] = readout("Rgr") * ln2 * vslice[Z_] * vslice[P_] / (readout("K") + vslice[P_]);
    F[find("Pmort")] = m * vslice[P_];
    F[find("Zmort")] = readout("g") * vslice[Z_];
    F[find("excret")] = readout("Rex") * F[find("ingest")]; // note: this may not be how ROMS defines it
    
    return F; 
  }
  
  void adjustments_after_calc_step() {
    totalN = vars[N_].current + vars[P_].current + vars[Z_].current;
    for (int i=0; i<nvars; i++) {
      vars[i].min = 0.0001;
      vars[i].max = totalN;
      vars[i].scale = totalN0/2;
    }
  }
  
} // NPZsystem


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


class NNPZDsystem extends dynamicalSystem {
  
  float totalN, totalN0; // internal variables
  int NO3_, NH4_, P_, Z_, D_, e0_, e1_;
  
  NNPZDsystem() { // separating the constructor and define() means you can easily make a variant that calls the same define() but does other things in the constructor
    define();
  }  
  
  dynamicalSystem clone() {
    return new NNPZDsystem();
  }
  
  void define() {  
    int numFramesToSave = 100;
    int maxvars = 7;
    int maxfluxes = 12;
    int maxparams = 14;
    allocate(numFramesToSave, maxvars, maxfluxes, maxparams);
    name = "NNPZD";
    shortname = "NNPZD";
    totalN0 = 35;
    float P0 = 1; float Z0 = 1; float D0 = 0;
    modelTimePerSecond = 2;
    timeUnits = "days";
    varUnits = "µmol/L";
    tolerance = 1e-5;

    color Pcol = greenColor; color Zcol = redColor;
    color Dcol = color(60,70,130); color NH4col = color(20,30,90); color NO3col = color(0,0,50);
    float R = 1; float po = PI/180;
    NO3_ = addVar("nitrate","NO3", totalN0-P0-Z0-D0, NO3col, sqrt(3)*R*cos(225*po)/2,(sqrt(3)*R*sin(225*po)-R)/2);
    NH4_ = addVar("ammonium","NH4", 0, NH4col, 0, -R);
    P_ = addVar("phytoplankton","P", P0, Pcol,R*cos(150*po),R*sin(150*po));
    Z_ = addVar("zooplankton","Z", Z0, Zcol,R*cos(30*po),R*sin(30*po));
    D_ = addVar("detritus","D",D0, Dcol,R*cos(-30*po),R*sin(-30*po));
    e0_ = addEndpoint("NO3 source","e0", vars[NO3_].xc,vars[NO3_].yc-0.5*R);
    e1_ = addEndpoint("D sink","e1", vars[D_].xc,vars[D_].yc-0.5*R);
    
    addFlux("NO3 uptake","NO3upt", NO3_, P_, Pcol,120-20);
    addFlux("NH4 uptake","NH4upt", NH4_, P_, Pcol,120);
    addFlux("ingestion","ingest", P_, Z_, Zcol,120);
    addFlux("egestion","egest", P_, D_, Zcol,180-30);
    addFlux("P mortality to D","PmortD", P_, D_, Dcol,180-40-10);
    addFlux("P mortality to N","PmortN", P_, NH4_, NH4col,120);
    addFlux("Z mortality","Zmort", Z_, D_, Dcol,60);
    addFlux("ingestion-related excretion","excret", Z_, NH4_, Zcol,120-30);
    addFlux("basal metabolism","metab", Z_, NH4_, NH4col,120-30-10);
    addFlux("remineralization","remin", D_, NH4_, NH4col,60);
    addFlux("NO3 supply","supply", e0_,NO3_, NO3col,0);
    addFlux("detritus sinking","sink", D_,e1_, Dcol,0);
    
    addParam("NO3 supply (µmol/L/day)","S",1, 0,1, NO3col);
    addParam("P max growth (1/day)","Vm",1.5,0,2, Pcol);
    addParam("Ks NO3 (µmol/L)","Kno3",1, 0, 1, Pcol);
    addParam("Ks NH4 (µmol/L)","Knh4",1, 0, 1, Pcol);
    addParam("P mortality (1/day)","m",0.1, 0, 0.4, Pcol);
    addParam("P mort. frac. to N","f",0.5, 0, 1, Pcol);
    addParam("Z max growth (1/day)","Rgr",0.6, 0, 2, Zcol);
    addParam("Ks grazing (µmol/L)","Kphy",2, 0, 4, Zcol);
    addParam("egestion fraction","gamma",0.3, 0, 1, Zcol);
    addParam("excretion fraction","Rex",0.2, 0, 1, Zcol);
    addParam("Z mortality (1/day)","g",0.1, 0, 0.4, Dcol);
    addParam("Z basal metabolism (1/day)","Rbm",0.1, 0, 0.4, Dcol);
    addParam("D remineralization (1/day)","r",0.1, 0, 1, Dcol);
    addParam("D sinking (1/day)","rs",0.3, 0, 0.4, Dcol);
  }
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    float Kno3 = readout("Kno3"); float Knh4 = readout("Knh4"); float Vm = readout("Vm"); float gamma = readout("gamma");
    float f = readout("f"); float m = readout("m");
    
    F[find("NO3upt")] = Vm * vslice[P_] * vslice[NO3_] / (Kno3 + vslice[NO3_] + Kno3/Knh4*vslice[NH4_]);
    F[find("NH4upt")] = Vm * vslice[P_] * Kno3/Knh4*vslice[NH4_] / (Kno3 + vslice[NO3_] + Kno3/Knh4*vslice[NH4_]);
    F[find("ingest")] = (1-gamma) * readout("Rgr") * vslice[Z_] * vslice[P_] / (readout("Kphy") + vslice[P_]);
    F[find("egest")] = gamma / (1-gamma) * F[find("ingest")];
    F[find("PmortN")] = f * m * vslice[P_];
    F[find("PmortD")] = (1-f) * m * vslice[P_];
    F[find("Zmort")] = readout("g") * vslice[Z_];
    F[find("remin")] = readout("r") * vslice[D_];
    F[find("excret")] = readout("Rex") * F[find("ingest")]; // note: this may not be how ROMS defines it
    F[find("metab")] = readout("Rbm") * vslice[Z_];
    F[find("supply")] = readout("S");
    F[find("sink")] = readout("rs") * vslice[D_];
    
    return F; 
  }
  
  void adjustments_after_calc_step() {
    totalN = vars[NO3_].current + vars[NH4_].current + vars[P_].current + vars[Z_].current + vars[D_].current;
    for (int i=0; i<nvars; i++) {
      vars[i].min = 0.0001;
      vars[i].max = totalN;
      vars[i].scale = totalN0/2;
    }
  }
  
} // NNPZDsystem




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




class allo_upwelling_system1D extends dynamicalSystem {

  int NP, NZ, Nsub; // number of phytos, zooplankton, and subsystems
  int iN0; int[] iT, iN, iD; int[][] iP, iZ; // statevar indices
  int iFsupp; int[] iFremin, iFsink; int[][] iFup, iFsen, iFeg, iFresp, iFadv; int[][][] iFing; // flux indices
  float[] x, mu, sigma, k, K, rho, thresh, Imax, rpp; float[][] pref, scaledstock; // internal param arrays
  
  float totalN0, totalN;
  
  float vang, sinvang, cosvang; // screen layout parameters
  


  allo_upwelling_system1D(int nphytos, int nzoop, int nsubsystems, float totalnitrogen) {

    name = "allometric/upwelling";
    shortname = "allo-upwelling";
    modelTimePerSecond = 2;
    timeUnits = "days";
    varUnits = "µM";
    tolerance = 1e-5;
    int numFramesToSave = 100;
    int params_to_allocate = 50;
    if (nsubsystems*(nphytos+nzoop)>40) {showFluxes = false;}
    showHistory = false;
    showDetails = false;
    showDetailsButton = true;

    // set values passed to the constructor;
    // allocate arrays
    totalN0 = totalnitrogen;
    totalN = totalN0;
    NP = nphytos;
    NZ = nzoop;
    Nsub = nsubsystems;
    int vars_to_allocate = Nsub*(NP+NZ+3) + 1 +10;
    int fluxes_to_allocate = Nsub*(NP+NZ)*(NZ+4)+10+100;
    allocate(numFramesToSave, vars_to_allocate, fluxes_to_allocate, params_to_allocate);

    // layout parameters
    dot_minsize = dot_minsize/2;
    dot_stdsize = 0.25-0.03*(Nsub-2);
    vang = 15;
    float spacing = 0.7+0.14*(Nsub-2);
    color[] Pcolor = new color[NP];
    color[] Zcolor = new color[NZ];
    for (int i=0; i<NP; i++) {Pcolor[i] = randomshift(greenColor);}
    for (int i=0; i<NZ; i++) {Zcolor[i] = randomshift(redColor);}
    color Dcolor = lighten(blueColor);
    color Ncolor = blueColor;
    color N0color = darken(blueColor);
        
    // create & lay out subsystems
    float ystart = 0;
    float xstart = 0;
    if (Nsub==1) {
      ystart=-0.1;
      xstart=-0.7;
    }
    cosvang = cos(vang/180*PI); sinvang = sin(vang/180*PI);
    iN0 = addVar("deep NO3","deep N", totalN0, N0color, xstart+spacing-sinvang, ystart-cosvang); // create deep N pool
    iT = new int[Nsub];
    float dx = 2.0*spacing/(float)(constrain(Nsub,2,20)-1);
    for (int ns=0; ns<Nsub; ns++) {
      iT[ns] = addVar("total N subsystem "+str(ns+1),"T"+str(ns+1), 0, dkGrayColor, xstart+spacing-ns*dx, ystart); // create hidden variables for totals in each subsystem
      vars[iT[ns]].display = false;
      vars[iT[ns]].constant = true;
    }
    
    // create other variables
    // note that positions are all assigned dynamically, and are nan here
    iP = new int[Nsub][NP];
    iZ = new int[Nsub][NZ];
    iN = new int[Nsub];
    iD = new int[Nsub];
    for (int ns=0; ns<Nsub; ns++) {
      iN[ns] = addVar("nutrients "+str(ns+1), "N", 0, Ncolor, nan, nan);
      iD[ns] = addVar("detritus "+str(ns+1), "D", 0, Dcolor, nan, nan); 
      for (int i=0; i<NP; i++) {
        iP[ns][i] = addVar("phyto "+str(ns+1)+","+str(i), "", 0.01, Pcolor[i], nan, nan);
      }
      for (int i=0; i<NZ; i++) {
        iZ[ns][i] = addVar("zoop "+str(ns+1)+","+str(i), "", 0.01, Zcolor[i], nan, nan);
      }
    }
    
    // create fluxes
    
    iFup = new int[Nsub][NP];
    iFsen = new int[Nsub][NP];
    iFeg = new int[Nsub][NP+NZ];
    iFresp = new int[Nsub][NZ];
    iFing = new int[Nsub][NP+NZ][NZ];
    iFadv = new int[Nsub-1][NP+NZ+2];
    iFremin = new int[Nsub];
    iFsink = new int[Nsub];
    
    for (int ns=0; ns<Nsub; ns++) {
      // phyto solo activities
      for (int i=0; i<NP; i++) {
        String indices = str(ns+1)+","+str(i);
        iFup[ns][i] = addFlux("uptake "+indices, "up"+indices, iN[ns], iP[ns][i], Pcolor[i], 90);
        iFsen[ns][i] = addFlux("senescence "+indices, "sen"+indices, iP[ns][i], iD[ns], Dcolor, -20);
        iFeg[ns][i] = addFlux("egestionP "+indices, "egp"+indices, iP[ns][i], iD[ns], Dcolor, 20); 
      }
      // zoop solo activities
      for (int i=0; i<NZ; i++) {
        String indices = str(ns+1)+","+str(i);
        iFresp[ns][i] = addFlux("respiration "+indices, "resp"+indices, iZ[ns][i], iN[ns], Dcolor, 20+i);
        iFeg[ns][NP+i] = addFlux("egestionZ "+indices, "egz"+indices, iZ[ns][i], iD[ns], Dcolor, 90);
      }
      // ingestion fluxes
      for (int j=0; j<NZ; j++) {
        for (int i=0; i<NP; i++) {
          String indices = str(ns+1)+","+str(i);
          iFing[ns][i][j] = addFlux("ingestionP"+indices, "ingp"+indices, iP[ns][i], iZ[ns][j], Zcolor[j], 90);
        }
        for (int i=0; i<NZ; i++) {
          String indices = str(ns+1)+","+str(i);
          iFing[ns][NP+i][j] = addFlux("ingestionZ"+indices, "ingz"+indices, iZ[ns][i], iZ[ns][j], Zcolor[j], -150);
        }
      }
      // remineralization & sinking
      iFremin[ns] = addFlux("remineralization "+str(ns+1), "remin"+str(ns+1), iD[ns], iN[ns], Ncolor, 90);
      float curv = 60; if (ns>0) {curv = -60;}
      iFsink[ns] = addFlux("sinking "+str(ns+1), "sink"+str(ns+1), iD[ns], iN0, N0color, curv);
      // advection between subsystems
      if (ns<Nsub-1) {
        for (int i=0; i<NP; i++) {
          String indices = str(ns+1)+","+str(i);
          iFadv[ns][i] = addFlux("advectionP "+indices, "advp"+indices, iP[ns][i], iP[ns+1][i], grayColor, 1);
        }
        for (int i=0; i<NZ; i++) {
          String indices = str(ns+1)+","+str(i);
          iFadv[ns][NP+i] = addFlux("advectionZ "+indices, "advz"+indices, iZ[ns][i], iZ[ns+1][i], grayColor, 1);
        }
        iFadv[ns][NP+NZ+0] = addFlux("advectionN "+str(ns+1), "advn"+str(ns+1), iN[ns], iN[ns+1], grayColor, 1);
        iFadv[ns][NP+NZ+1] = addFlux("advectionD "+str(ns+1), "advd"+str(ns+1), iD[ns], iD[ns+1], grayColor, 1);        
      }
    }  
    // deep nutrient supply
    iFsupp = addFlux("supply 1", "supp1", iN0, iN[0], blueColor, -5);
    fluxes[iFsupp].arrowheadpos = 0.5;

    
    // add normal parameters
    if (Nsub>1) {
      addParam("surface box width (km)","b",20,5,100, darken(blueColor));
      addParam("cross-shelf velocity (cm/s)","u",5, 0,25, darken(blueColor));
    } else {
      addConstant("surface box width (km)","b",20);
      addConstant("cross-shelf velocity (cm/s)","u",5);
    } 
    addParam("N supply intensity","supp_intensity", 1, 0, 10, darken(blueColor));
    addParam("upwelling-downwelling frequency (cycles/d)",  "omega", 0, 0,0.15, darken(blueColor));
    addParam("D sinking rate (1/d)", "woh",8.0/20, 0,1, darken(blueColor)); // w_sink / h
    addParam("D remineralization rate (1/d)","rem",0.1, 0,0.3, blueColor);
    addConstant("min phyto size (µm)",     "xpmin",1, 1,20, greenColor);
    addConstant("max phyto size (µm)",     "xpmax",10, 1,40, greenColor);
    addConstant("min zoop size (µm)",     "xzmin",3, 3,150, redColor);
    addConstant("max zoop size (µm)",     "xzmax",2000, 3, 2000, redColor);
    addConstant("pred-prey size ratio (1 µm)", "rpp0", 2);
    addConstant("pred-prey size ratio exponent", "expo_rpp", 0.4); 
    addParam("spread in pred-prey size ratio","dlogrpp",0.23,  0.0,0.3, redColor); // preference falls off 1/e this far from rppmean
    addConstant("P metabolic exponent",      "expo_metab_P",-0.45);
    addConstant("P max growth (1 µm)",     "mu0",1.95);
    addConstant("N->P half-sat (1 µm)",    "k0",0.08);
    addConstant("N->P half-sat exponent",  "expo_k",1);
    addConstant("Z metabolic exponent",      "expo_metab_Z",-0.39);
    addConstant("Z max ingest (1 µm)",     "Imax0",18.5);
    addParam("PZ->Z half-sat (1 µm)",   "K0",3,0.6,15,redColor);
    addConstant("PZ->Z half-sat exponent", "expo_K",0);
    addConstant("Z respiration (1 µm)",    "rho0",1.85);
    addParam("P mortality (fractional)",      "sigma",0.1,0,1,blueColor);
    addParam("egestion (fractional)",          "gamma",0.3,  0,1, redColor);
    addParam("relative grazing threshhold","theta",0.2,0,0.3,redColor);
    addConstant("spread in mu",            "dlogmu",0);
    addConstant("spread in k",             "dlogk",0);
    addConstant("spread in Imax",          "dlogImax",0);
    addConstant("spread in K",             "dlogK",0);
    addConstant("spread in rho",           "dlogrho",0);
    addConstant("spread in theta",         "dlogtheta",0);
    addConstant("spread in gamma",         "dloggamma",0);  
    addParam("Z preference for Z",         "zeta",1, 0,1, redColor);

    // allocate internal parameters
    x = new float[NP+NZ]; // P,Z size
    mu = new float[NP]; // phyto max growth
    sigma = new float[NP]; // senescence as fraction of mu
    k = new float[NP]; // half-sat for phyto uptake
    K = new float[NP+NZ]; // half-sat when being eaten
    rho = new float[NZ]; // Z respiration rate
    thresh = new float[NZ]; // Z grazing threshhold
    Imax = new float[NZ]; // Z max ingestion
    rpp = new float[NZ]; // optimal predator-prey size ratio
    pref = new float[NP+NZ][NZ]; // feeding preferences
    scaledstock = new float[NP+NZ][NZ]; // scratch variable for calcfluxes
   
    RESPOND_TO_PARAMCHANGE();
    
    adjustments_after_calc_step();
    float extra = totalN-totalN0;
    vars[iN0].current = vars[iN0].current - extra;
    vars[iN0].initial = vars[iN0].initial - extra;
    
    // add a few more variables to serve as a scale bar
    int scaleval = 5;
    if (Nsub==1) {scaleval=1;}
    int q = addVar("scalebar "+str(scaleval),str(scaleval)+" mmol/m3",scaleval, grayColor, -0.2, -1.5);
    vars[q].constant = true;
    
       
  } // define allo_upwelling
  
  
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    
    float sink_rate = readout("woh");
    float remin_rate = readout("rem");
    float zeta = readout("zeta");
    float gamma = readout("gamma");
    float go1mg = gamma/(1-gamma);
//  float uob = readout("uob");
    float uob = readout("u")*100*1000/86400/readout("b");
//  float uob_now = uob * sq(cos(3.14159*readout("omega")*tcurrent)); // upwelling - relaxation
    float uob_now = uob * constrain(cos(2*3.14159*readout("omega")*tcurrent),0,big);  // upwelling - downwelling

    // N supply
    F[iFsupp] = uob_now * readout("supp_intensity") * vslice[iN0];
    for (int ns=0; ns<Nsub; ns++) { 
      // uptake & senescence
      for (int i=0; i<NP; i++) {
        F[iFup[ns][i]] = mu[i] * vslice[iN[ns]]/(k[i]+vslice[iN[ns]]) * vslice[iP[ns][i]];
        F[iFsen[ns][i]] = sigma[i] * mu[i] * vslice[iP[ns][i]];
      }   
      // respiration
      for (int i=0; i<NZ; i++) {
        F[iFresp[ns][i]] = rho[i] * vslice[iZ[ns][i]];
      }
      // ingestion
      float totalfood, totalingest, frac; 
      for (int j=0; j<NZ; j++) {
        totalfood = 0;
        for (int i=0; i<NP; i++) {
          scaledstock[i][j] = pref[i][j] * vslice[iP[ns][i]] / K[i];
          totalfood = totalfood + scaledstock[i][j];
        }
        for (int i=0; i<NZ; i++) {
          scaledstock[NP+i][j] = zeta * pref[NP+i][j] * vslice[iZ[ns][i]] / K[NP+i];
          totalfood = totalfood + scaledstock[NP+i][j];
        }
        // now that we know the totalfood for grazer j...
        totalingest = Imax[j] * vslice[iZ[ns][j]] * constrain(totalfood - thresh[j], 0,big) / (1 + totalfood - thresh[j]);
        for (int i=0; i<NP+NZ; i++) {
          if (totalfood > thresh[j]) {
            frac = constrain(scaledstock[i][j] / (totalfood-thresh[j]), 0,1);
            F[iFing[ns][i][j]] = frac * totalingest;
          } else {
            F[iFing[ns][i][j]] = 0;
          }
        }
      } // end loop over Z
      // remineralization
      F[iFremin[ns]] = remin_rate * vslice[iD[ns]];
      F[iFsink[ns]] = sink_rate * vslice[iD[ns]];
      // egestion: for each prey i, add up contributions from all grazers j
      for (int i=0; i<NP+NZ; i++) {
        F[iFeg[ns][i]] = 0;
        for (int j=0; j<NZ; j++) {
          F[iFeg[ns][i]] = F[iFeg[ns][i]] + go1mg * F[iFing[ns][i][j]];
        }
      }
      // advection
     if (ns<Nsub-1) {
        for (int i=0; i<NP; i++) {F[iFadv[ns][i]] = uob_now * vslice[iP[ns][i]];}
        for (int i=0; i<NZ; i++) {F[iFadv[ns][NP+i]] = uob_now * vslice[iZ[ns][i]];}
        F[iFadv[ns][NP+NZ+0]] = uob_now * vslice[iN[ns]];
        F[iFadv[ns][NP+NZ+1]] = uob_now * vslice[iD[ns]];
      }
    } // end loop over subsystems
                
    return F; 
  } // calcfluxes for allo-upwelling
  
  
  
  void adjustments_after_calc_step() {
    totalN = vars[iN0].current;
    for (int ns=0; ns<Nsub; ns++) {
      for (int i=0; i<NP; i++) {totalN = totalN + vars[iP[ns][i]].current;}
      for (int i=0; i<NZ; i++) {totalN = totalN + vars[iZ[ns][i]].current;}
      totalN = totalN + vars[iN[ns]].current;
      totalN = totalN + vars[iD[ns]].current;
    }
    for (int i=0; i<nvars; i++) {
      vars[i].min = 0.0001;
      vars[i].max = totalN;
      vars[i].scale = totalN0/2;
      vars[i].current = constrain(vars[i].current, vars[i].min, vars[i].max);
    }
  }
  
  
  
  void RESPOND_TO_PARAMCHANGE() {

    float xpmin = readout("xpmin");    float logxpmin = log(xpmin)/log(10);
    float xpmax = readout("xpmax");    float logxpmax = log(xpmax)/log(10);
    float xzmin = readout("xzmin");    float logxzmin = log(xzmin)/log(10);
    float xzmax = readout("xzmax");    float logxzmax = log(xzmax)/log(10);
    for (int i=0; i<NP; i++) {
      x[i] = pow(10,logxpmin + (logxpmax-logxpmin)*(float)i/((float)NP-1)); // set phyto sizes
      for (int ns=0; ns<Nsub; ns++) {vars[iP[ns][i]].shortname = str_sigfigs(x[i]);} // set phyto names
    }
    for (int i=0; i<NZ; i++) {
      x[NP+i] = pow(10,logxzmin + (logxzmax-logxzmin)*(float)i/((float)NZ-1)); // set zoop sizes
      for (int ns=0; ns<Nsub; ns++) {vars[iZ[ns][i]].shortname = str_sigfigs(x[NP+i]);} // set zoop names
    }
    
    for (int i=0; i<NP; i++) { // phyto size-based params
      mu[i] = readout("mu0") * pow(x[i],readout("expo_metab_P"));
      sigma[i] = readout("sigma");
      k[i] = readout("k0") * pow(x[i],readout("expo_k"));
    }    
    for (int i=0; i<NZ; i++) { // zoop size-based params
      rho[i] = readout("rho0") * pow(x[NP+i],readout("expo_metab_Z"));
      Imax[i] = readout("Imax0") * pow(x[NP+i],readout("expo_metab_Z"));
      thresh[i] = readout("theta");
      rpp[i] = readout("rpp0") * pow(x[NP+i],readout("expo_rpp"));
    }    
    float dlogrpp = readout("dlogrpp"); // params that apply to all prey
    for (int i=0; i<NP+NZ; i++) {
      K[i] = readout("K0") * pow(x[i],readout("expo_K"));
      for (int j=0; j<NZ; j++) {
        pref[i][j] = exp(-sq(log(x[NP+j]/x[i]/rpp[j])/log(10)/dlogrpp)); // preference of j for i, as function of size difference
        pref[i][j] = constrain(pref[i][j],0,1);
      }
    }
  }
  
  
  
  void adjustments_before_drawing() { // set dot positions
    float x0,y0,R;
    for (int ns=0; ns<Nsub; ns++) {
      x0 = vars[iT[ns]].xc;
      y0 = vars[iT[ns]].yc;
      
      R = dot_stdsize + vars[iN[ns]].radius;
      vars[iN[ns]].xc = x0 - R * cosvang;
      vars[iN[ns]].yc = y0 - R * sinvang;
      
      R = dot_stdsize + vars[iD[ns]].radius;
      vars[iD[ns]].xc = x0 + R * sinvang;
      vars[iD[ns]].yc = y0 - R * cosvang;
      
      R = dot_stdsize;
      for (int i=0; i<NP; i++) {
        R = R + vars[iP[ns][i]].radius;
        vars[iP[ns][i]].xc = x0 - R * sinvang;
        vars[iP[ns][i]].yc = y0 + R * cosvang;
        R = R + vars[iP[ns][i]].radius;        
      }
      
      R = dot_stdsize;
      for (int i=0; i<NZ; i++) {
        R = R + vars[iZ[ns][i]].radius;
        vars[iZ[ns][i]].xc = x0 + R * cosvang;
        vars[iZ[ns][i]].yc = y0 + R * sinvang;
        R = R + vars[iZ[ns][i]].radius;        
      }      
    }
    
    for (int i=0; i<nfluxes; i++) {
      // hide fluxes that = 0 instead of drawing them in white
      if (fluxes[i].current < tolerance*fluxes[i].scale) {fluxes[i].current = nan;}
      fluxes[i].width = fluxes[i].width/2;
    }
    
  }
  
  
  void drawonevarlabel(int i, float xcen, float ycen, float r) { // override so that small labels aren't drawn
    color col = vars[i].thecolor;
    if (vars[i].current < tolerance*vars[i].scale) {col = 255;}
    textAlign(CENTER);
    if (r>smallfontsize) {
      setFont("flexi");
      fill(lighten(lighten(lighten(col))));
      textSize(r);
      textVAlign(vars[i].shortname,xcen,ycen,r,"middle");
    } else {
//      if (col != 255) {fill(darken(col));} else {fill(255);}
//      setFont("small");
//      textVAlign(vars[i].shortname,xcen,ycen-r-1,smallfontsize,"bottom");
    }
  }
 
 
  void draw_other_things() {
    float xscale = networkRect[2]/2;
    float xoffset = networkRect[0]+xscale;
    float yscale = networkRect[3]/2;
    float yoffset = networkRect[1]+yscale;
    if (showDetails) {
      label_boxes(xscale, yscale, xoffset, yoffset);
    }
    if (showDetails) {
      label_statevars(xscale, yscale, xoffset, yoffset);
    }       
    if (showDetails) {
      displayPreferences();
    }
  }
  
  
  void label_boxes(float xscale, float yscale, float xoffset, float yoffset) {
    float x,y;
    String inshoreStr;
    setFont("big");
    fill(dkGrayColor);
    textAlign(LEFT);
    x = vars[iN0].xc + 1.3 * vars[iN0].radius;
    y = vars[iN0].yc;
    textVAlign("subsurface nutrient pool",x*xscale+xoffset+30,-y*yscale+yoffset,bigfontsize,"middle"); // the +30 is a kluge because "deep N" is a wide label
    if (Nsub>1) {inshoreStr = "surface layer, nearshore end";} else {inshoreStr = "surface layer";}
    textAlign(LEFT);
    x = vars[iD[0]].xc + 0.9 * vars[iD[0]].radius;
    y = vars[iD[0]].yc - 0.9 * vars[iD[0]].radius;
    textVAlign(inshoreStr,x*xscale+xoffset,-y*yscale+yoffset+bigfontsize,bigfontsize,"top");
    if (Nsub>1) {
      textAlign(RIGHT);
      x = vars[iD[Nsub-1]].xc - vars[iD[Nsub-1]].radius;
      y = vars[iD[Nsub-1]].yc - vars[iD[Nsub-1]].radius;
      textVAlign("surface layer, offshore end",x*xscale+xoffset,-y*yscale+yoffset,bigfontsize,"top");
      
      float dx = (vars[iT[0]].xc - vars[iT[1]].xc)/2.0 * xscale;
      x = vars[iT[1]].xc * xscale + xoffset;
      y = -vars[iP[1][NP-1]].yc * yscale + yoffset - bigfontsize;
      stroke(dkGrayColor);
      strokeWeight(2);
      line(x-dx,y,x+dx,y);
      line(x-dx,y-4,x-dx,y+4);
      line(x+dx,y-4,x+dx,y+4);
      noStroke();
      textAlign(CENTER); textVAlign(str(round(readout("b"))) + " km",x,y-4,bigfontsize,"bottom");
    }
    
  }
  
  void label_statevars(float xscale, float yscale, float xoffset, float yoffset) {
    float x,y;
    setFont("med");
    fill(greenColor);
    x = vars[iP[0][NP-1]].xc;
    y = vars[iP[0][NP-1]].yc + vars[iP[0][NP-1]].radius;
    textAlign(CENTER);
    textVAlign("phytoplankton (organism size in µm)",x*xscale+xoffset,-y*yscale+yoffset,medfontsize,"bottom");
    fill(redColor);
    x = vars[iZ[0][0]].xc + 1.2*vars[iZ[0][0]].radius;
    y = vars[iZ[0][0]].yc - 1.2*vars[iZ[0][0]].radius;
    textAlign(LEFT);
    textVAlign("zooplankton (organism size in µm)",x*xscale+xoffset,-y*yscale+yoffset,medfontsize,"top");
    fill(vars[iN[0]].thecolor);
    x = vars[iN[0]].xc - 0.8 * vars[iN[0]].radius;
    y = vars[iN[0]].yc - 0.8 * vars[iN[0]].radius;
    textAlign(RIGHT);
    textVAlign("available nutrients",x*xscale+xoffset,-y*yscale+yoffset,medfontsize,"middle");
    fill(vars[iD[0]].thecolor);
    x = vars[iD[0]].xc + 0.9 * vars[iD[0]].radius;
    y = vars[iD[0]].yc - 0.9 * vars[iD[0]].radius;
    textAlign(LEFT);
    textVAlign("detritus",x*xscale+xoffset,-y*yscale+yoffset,medfontsize,"middle");
  }
  
  
  void displayPreferences() {
    float[][] pr = new float[NZ][NP+NZ];
    float zeta = readout("zeta");
    for (int i=0; i<NP+NZ; i++) {
      for (int j=0; j<NZ; j++) {
        if (i>=NP) {
          pr[j][i] = log(zeta*pref[i][j]);
        } else {
          pr[j][i] = log(pref[i][j]);
        }
      }
    }
    int dy = round(historyRect[3]*0.67);
    int dx = round((NP+NZ)/(float)NZ*dy);
    if (dx > historyRect[2]) {
      dx = historyRect[2];
      dy = round(NZ/((float)NP+NZ)*dx);
    }
    int x0 = windowWidth - historyRect[0] - historyRect[2];
    int y0 = historyRect[1] + round(historyRect[3]*0.33);
    colormappedgrid(pr, log(tolerance), 0, lighten(ltGrayColor), dkGrayColor, defineRect(x0,y0,dx,dy));
    noStroke();
    fill(dkGrayColor);
    setFont("med");
    textAlign(CENTER);
    textVAlign("grazing preferences",x0+dx/2, y0-6, medfontsize,"bottom");
/*    textVAlign("P",x0+NP/((float)NP+NZ)/2*dx, y0+dy+5,medfontsize,"top");
    textVAlign("Z",x0+dx-NZ/((float)NP+NZ)/2*dx, y0+dy+5, medfontsize,"top");
    textAlign(RIGHT);
    textVAlign("Z ",x0,y0+dy/2,medfontsize,"middle");  */
    textAlign(LEFT); textVAlign("P size >>",x0, y0+dy+5,medfontsize,"top");
    textAlign(LEFT); textVAlign("Z size >>",x0+dx*NP/((float)NP+NZ), y0+dy+5, medfontsize,"top");
    textAlign(RIGHT); textVAlign("Z size ",x0, y0+dy/2, medfontsize,"middle");
  }
    
     
}



class romsBio2p2z extends dynamicalSystem {
  
  float totalN, totalN0; // internal variables
  int NO3_, NH4_, P_, P2_, Z_, Z2_, D_, par_, scale_;
  
  romsBio2p2z() {
    define();
  }  
  
  dynamicalSystem clone() {
    return new romsBio2p2z();
  }
  
  void define() {  
    adjustFluxesAtMaxMin = true;
    int numFramesToSave = 100;
    int maxvars = 11;
    int maxfluxes = 17;
    int maxparams = 40;
    allocate(numFramesToSave, maxvars, maxfluxes, maxparams);
    name = "roms PPZZ";
    shortname = "2p2z";
    modelTimePerSecond = 2;
    timeUnits = "days";
    varUnits = "µmol/L";
    tolerance = 1e-5;
    dot_stdsize *= 3;

    color Pcol = greenColor;
    color P2col = darken(Pcol);
    color Zcol = redColor;
    color Z2col = darken(Zcol);
    color Dcol = color(60,70,130);
    color NH4col = color(20,30,90);
    color NO3col = color(0,0,50);
    
    float sp = 0.7;
    NO3_ = addVar("nitrate","NO3", 0, NO3col, -sp, -sp);
    NH4_ = addVar("ammonium","NH4", 0, NH4col, 0, -sp);
    D_ = addVar("detritus","D",0, Dcol, sp, -sp);
    P_ = addVar("phytoplankton","P", 0, Pcol, -sp/2, 0);
    P2_ = addVar("phytoplankton2","P2", 0, P2col, -sp*0.866, sp);
    Z_ = addVar("zooplankton","Z", 0, Zcol, sp/2, 0);
    Z2_ = addVar("zooplankton2","Z2", 0, Z2col, sp*0.866, sp);
    par_ = addVar("PAR","PAR", 100, yellowColor, 0,0); // this is where communicationAfterCalcStep() stores PAR for each cell
    vars[par_].display = false;
    
    // scale bar: off by default
    color scaleColor = transparency(color(255),0.5);
    scale_ = addVar("scale","scale", 10, color(255), 2,-4);
    vars[scale_].display = false;
    vars[addVar("scale1","scale1",1,color(255),-2,-4)].display = false;
    vars[addVar("scale2","scale2",5,color(255),-0.2,-4)].display = false;
    
    setInitialConditions();
        
    addFlux("NO3 uptake","NO3upt", NO3_, P_, Pcol,60);
    addFlux("NO3 uptake 2","NO3upt2", NO3_, P2_, P2col,60);
    addFlux("NH4 uptake","NH4upt", NH4_, P_, Pcol,120);
    addFlux("NH4 uptake 2","NH4upt2", NH4_, P2_, P2col,120);
    addFlux("P mortality","Pmort", P_, D_, Dcol, -10);
    addFlux("P mortality2","Pmort2", P2_, D_, Dcol, -10);
    addFlux("ingestion","ingest", P_, Z_, Zcol,120);
    addFlux("ingestion2","ingest2", P2_, Z2_, Z2col,110);
    addFlux("egestion","egest", P_, D_, Zcol,10);
    addFlux("egestion2","egest2", P2_, D_, Z2col,10);
    addFlux("ingestion-related excretion","excret", Z_, NH4_, Zcol,120);
    addFlux("ingestion-related excretion2","excret2", Z2_, NH4_, Z2col,120);
    addFlux("basal metabolism","metab", Z_, NH4_, NH4col,80);
    addFlux("basal metabolism2","metab2", Z2_, NH4_, NH4col,80);
    addFlux("Z mortality","Zmort", Z_, D_, Dcol,60);
    addFlux("Z mortality2","Zmort2", Z2_, D_, Dcol,60);
    addFlux("remineralization","remin", D_, NH4_, NH4col,60);

    addConstant("attenuation by seawater","attSW", 0.15);
    addConstant("attenuation by chl","attChl", 0.013);
    addConstant("PAR fraction","parFrac", 0.43);
    addConstant("chl:N ratio","chl2N", 1); // this is handled dynamically in roms
    addConstant("light-growth initial slope","alpha",0.025); // in roms, PhyIS
    addParam("P max growth (1/day)","Vm",1.4, 0,2, Pcol);
    addParam("P2 max growth (1/day)","Vm2",1.4, 0,2, P2col);
    addConstant("Ks for uptake (µmol/L)","k",0.08, 0.08, 2, Pcol); // the inverse of this appears in bioFasham.h
    addConstant("Ks for uptake 2 (µmol/L)","k2",0.8, 0.08, 2, P2col);
    addParam("P, P2 mortality (1/day)","mp",0.05, 0, 0.4, Pcol);
    addParam("Z max grazing (1/day)","Imax",10, 0,10, Zcol);
    addParam("Z2 max grazing (1/day)","Imax2",1.2, 0, 4, Z2col);
    addParam("Ks for grazing (µmol/L)","K",3, 0.1, 10, Zcol); // in roms, the SQUARE of this is specified (...P^2/(K+P^2)...)
    addParam("egestion fraction","gamma",0.3, 0, 1, Zcol); // in roms, 1-ZooAE_N
    addParam("excretion fraction ","excr",0.4, 0, 1, Zcol);
    addConstant("Z basal metabolism (1/day)","basal",0);
    addConstant("Z2 basal metabolism (1/day)","basal2",0);
    addParam("Z,Z2 mortality (1/day)","mz",0.3, 0, 1.2, Zcol);
    addParam("D remineralization (1/day)","rem",0.1, 0, 1, Dcol);
    addParam("D sinking (m/day)","wsinkD",8, 0, 8, Dcol);
    addParam("P2 sinking (m/day)","wsinkP2",0, 0, 8, P2col);
    
    RESPOND_TO_PARAMCHANGE();
    
  }
  
  void setInitialConditions() {setInitialConditions(vars[scale_].current);}
  void setInitialConditions(float romsNO3) {
    vars[P_].initial = 0.1;
    vars[P2_].initial = 0.1;
    vars[Z_].initial = 0.01;
    vars[Z2_].initial = 0.01;
    vars[NO3_].initial = romsNO3;
    vars[NH4_].initial = 0;
    vars[D_].initial = 0;
  }
  
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    
    float P = vslice[P_];
    float P2 = vslice[P2_];
    float Z = vslice[Z_];
    float Z2 = vslice[Z2_];
    float NO3 = vslice[NO3_];
    float NH4 = vslice[NH4_];
    float D = vslice[D_];
    float PAR = vslice[par_];
    
    // 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("NO3upt")] = mu * NO3 / (k + NO3 + NH4) * P; // simplified nut-limitation expressions; assuming k_NO3 = k_NH4
    F[find("NH4upt")] = mu * NH4 / (k + NO3 + NH4) * P;
    float Vm2 = readout("Vm2");
    float mu2 = Vm2 * alphaI / sqrt(sq(Vm2) + sq(alphaI));
    float k2 = readout("k2") + tolerance;
    F[find("NO3upt2")] = mu2 * NO3 / (k2 + NO3 + NH4) * P2;
    F[find("NH4upt2")] = mu2 * NH4 / (k2 + NO3 + NH4) * P2;
    
    // 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;
      grazing2 = readout("Imax2") * P2 / (K + P2) * Z2;
    } else if (grazingResponse.equals("quadratic")) {
      grazing = readout("Imax") * sq(P) / (sq(K) + sq(P)) * Z;
      grazing2 = readout("Imax2") * sq(P2) / (sq(K) + sq(P2)) * Z2;
    }
    float gamma = readout("gamma");
    float excrFrac = readout("excr");
    F[find("ingest")] = (1-gamma) * grazing;
    F[find("ingest2")] = (1-gamma) * grazing2;
    F[find("egest")] = gamma * grazing;
    F[find("egest2")] = gamma * grazing2;
    F[find("excret")] = excrFrac * grazing;
    F[find("excret2")] = excrFrac * grazing2;
    
    // other loss terms
    F[find("metab")] = readout("basal") * Z;
    F[find("metab2")] = readout("basal2") * Z2;
    float mp = readout("mp");
    F[find("Pmort")] = mp * P;
    F[find("Pmort2")] = mp * P2;
    float mz = readout("mz");
    F[find("Zmort")] = mz * sq(Z);
    F[find("Zmort2")] = mz * sq(Z2);
    F[find("remin")] = readout("rem") * D;
    
    return F; 
  }
  
  void adjustments_after_calc_step() {
    totalN = vars[NO3_].current + vars[NH4_].current + vars[P_].current + vars[P2_].current + vars[Z_].current + vars[Z2_].current + vars[D_].current;
    for (int i=0; i<nvars; i++) {
      if (i != par_) {
        vars[i].min = 0.001;
        vars[i].max = totalN*100;
        vars[i].scale = vars[scale_].current;
      }
    }
  }
  
  void RESPOND_TO_PARAMCHANGE() {
    vars[D_].sinkingRate = readout("wsinkD");
    vars[P2_].sinkingRate = readout("wsinkP2");
  }
  
  float totalP() {
    return vars[P_].current + vars[P2_].current;
  }
  
  float totalNitrogen() {
    return totalN;
  }
  
  void showDiagnostics() {
    float Fup = fluxes[find("NO3upt")].current + fluxes[find("NH4upt")].current + fluxes[find("NO3upt2")].current + fluxes[find("NH4upt2")].current;
    float mu = Fup / totalP();
    float Fgraz = fluxes[find("ingest")].current + fluxes[find("egest")].current + fluxes[find("ingest2")].current + fluxes[find("egest2")].current;
    float g = Fgraz / totalP();
    textSize(40);
    textLeading(45);
    textAlign(CENTER);
    String diagnosticsString = "µ = " + str_sigfigs(mu) + "\n" + "g = " + str_sigfigs(g);
    fill(transparency(color(255),0.6));
    text(diagnosticsString,networkRect[0]+networkRect[2]/2,networkRect[1]+networkRect[3]/2);
    fill(transparency(yellowColor,0.6));
    text(diagnosticsString,networkRect[0]+networkRect[2]/2,networkRect[1]+networkRect[3]/2);
  }
  
} // romsBio2p2z




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



/* class NPDsystem extends dynamicalSystem {
  
  float totalN, totalN0; // internal variables
  int N_, P_, D_, par_, scale_;
  
  NPDsystem() {
    define();
  }  
  
  dynamicalSystem clone() {
    return new NPDsystem();
  }
  
  void define() {  
    adjustFluxesAtMaxMin = true;
    int numFramesToSave = 100;
    int maxvars = 11;
    int maxfluxes = 17;
    int maxparams = 40;
    allocate(numFramesToSave, maxvars, maxfluxes, maxparams);
    name = "NPD";
    shortname = "NPD";
    modelTimePerSecond = 2;
    timeUnits = "days";
    varUnits = "µmol/L";
    tolerance = 1e-5;
    dot_stdsize *= 3;

    color Pcol = greenColor;
    color Dcol = color(60,70,130);
    color Ncol = color(20,30,90);
    
    float sp = 0.7;
    N_ = addVar("nutrients","N", 0, Ncol, -sp, -sp);
    D_ = addVar("detritus","D",0, Dcol, sp, -sp);
    P_ = addVar("phytoplankton","P", 0, Pcol, 0, 0);
    par_ = addVar("PAR","PAR", 100, yellowColor, 0,0); // this is where communicationAfterCalcStep() stores PAR for each cell
    vars[par_].display = false;
    
    // scale bar: off by default
    color scaleColor = transparency(color(255),0.5);
    scale_ = addVar("scale","scale", 10, color(255), 2,-4);
    vars[scale_].display = false;
    vars[addVar("scale1","scale1",1,color(255),-2,-4)].display = false;
    vars[addVar("scale2","scale2",5,color(255),-0.2,-4)].display = false;
    
    setInitialConditions();
        
    addFlux("uptake","upt", N_, P_, Pcol,60);
    addFlux("P losses to N","lossN", P_, N_, Ncol, 60);
    addFlux("P losses to D","lossD", P_, D_, Dcol, 60);
    addFlux("remineralization","remin", D_, N_, Ncol,60);
    
    addConstant("attenuation by seawater","attSW", 0.15);
    addConstant("attenuation by chl","attChl", 0.013);
    addConstant("PAR fraction","parFrac", 0.43);
    addConstant("chl:N ratio","chl2N", 1); // this is handled dynamically in roms
    addConstant("light-growth initial slope","alpha",0.025); // in roms, PhyIS
    addParam("max uptake (1/day)","Vm",1.4, 0,2, Pcol);
    addParam("max grazing (1/day)","g",0.4, 0, 2, Pcol);
    addParam("Ks for uptake (µmol/L)","k",0.08, 0, 0.8, Pcol);
    addParam("Ks for grazing (µmol/L)","K",0, 0, 30, Dcol);
    addParam("mortality (1/day)","m",0.1, 0, 0.5, Pcol);
    addParam("P loss frac. to D", "Dfrac", 0.5, 0, 1, Dcol);
    addParam("D remineralization (1/day)","rem",0.1, 0, 1, Dcol);
    addParam("D sinking (m/day)","wsinkD",8, 0, 8, Dcol);
    
    RESPOND_TO_PARAMCHANGE();
    
  }
  
  void setInitialConditions() {setInitialConditions(vars[scale_].current);}
  void setInitialConditions(float romsNO3) {
    vars[P_].initial = 0.01;
    vars[N_].initial = romsNO3;
    vars[D_].initial = 0.01;
  }
  
  
  float[] CALCFLUXES(float[] vslice) {
    float[] F = new float[nfluxes];
    
    float P = vslice[P_];
    float N = vslice[N_];
    float D = vslice[D_];
    float PAR = vslice[par_];
    
    // 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;

    // phyto losses
    float K = readout("K") + tolerance;
    float losses = readout("g") * (P/(K+P)) * P;
    losses += readout("m") * P;
    float Dfrac = readout("Dfrac");
    F[find("lossN")] = (1-Dfrac) * losses;
    F[find("lossD")] = Dfrac * losses;
    
    // remineralization
    F[find("remin")] = readout("rem") * D;
    
    return F; 
  }
  
  void adjustments_after_calc_step() {
    totalN = vars[N_].current + vars[P_].current + vars[D_].current;
    for (int i=0; i<nvars; i++) {
      if (i != par_) {
        vars[i].min = 0.001;
        vars[i].max = totalN*100;
        vars[i].scale = vars[scale_].current;
      }
    }
  }
  
  void RESPOND_TO_PARAMCHANGE() {
    vars[D_].sinkingRate = readout("wsinkD");
  }
  
  float totalP() {
    return vars[P_].current;
  }
  
  float totalNitrogen() {
    return totalN;
  }
  
} // NPPDsystem
*/