class UpwellingConveyor extends Slice {

  UpwellingConveyor(Ecosystem eco, int M, int K, float Ltot, float Hcoast, float Hoc) {
    name = "2D upwelling conveyor";
    shortname = "upwelling";
    createDomain(M, K, Ltot, Hcoast, Hoc, 1, 1);
    linkToEcosystem(eco);
    
    internalTimestep = 0.001;

    Statevar mass = addVar("mass","mass","",Ncells);
    getVar("mass").setInitial(1);
    mass.advects = true;
    mass.mixes = true; 
    Statevar age = addVar("age","age","d",Ncells);
    age.advects = true;
    age.mixes = true;
    addFlux("aging","aging",null,age,Ncells);
    
    leftBC = new BoundaryCondition(this);
    leftBC.addVar("mass",1);    
    leftBC.addVar("age",0);
    
    allocateStorageForPhysics();

    addParam("surface PAR","PAR0","W/m2", 150, 0, 300);
    addParam("log vertical mixing","logkappa","m2/s", -4.5, -6, -4);
    addParam("D sinking rate", "wsinkD", "m/day", 10, 0, 100);
    addSinking("D", getParamVal("wsinkD"), 1);
    
    addParam("surface velocity", "u0", "m/s", -0.05, -0.3, 0.3);
    addParam("Ekman layer thickness", "hsurf", "m", 20, 0, 40);

    startArchive(1);
  }
  
  void updatePrivateParams() {
    float kappa0 = pow(10,getParamVal("logkappa")) * 86400;
    for (int m=0; m<M; m++) for (int k=0; k<K; k++) kappa[k][m] = kappa0;
    sinkingRate[getVarIndex("D")] = getParamVal("wsinkD");
    super.updatePrivateParams();
  }

  void calcFluxes() {
    for (int m=0; m<M; m++) {
      for (int k=0; k<K; k++) {
        getFlux("aging").current[km2ind(k,m)] = 1;
      }
    }
    super.calcFluxes();
  }
  
  void calcVelocity() {
    float u0 = getParamVal("u0") * 86400;
    float hsurf = getParamVal("hsurf");
    for (int m=0; m<M+1; m++) for (int k=0; k<K; k++) u[k][m] = 0;
    // apply velocity in the surface layer
    for (int m=0; m<M; m++) {
      for (int k=kbotu[m]; k<K; k++) {
        u[k][m] = u0 * exp(z_rho[k]/hsurf);
      }
    }
    // apply a compensating barotropic velocity
    for (int m=0; m<M+1; m++) {
      float uAtot = 0, Atot = 0;
      for (int k=kbotu[m]; k<K; k++) {
        uAtot += u[k][m] * Ayz[k][m];
        Atot += Ayz[k][m];
      }
      float umean = uAtot / Atot;
      for (int k=kbotu[m]; k<K; k++) u[k][m] -= umean;
    }
    super.calcVelocity();
  }


}