#include "all_vertical_richards.hpp"

namespace Kernel{
  void
  AllVerticalRichards::init(const Geometry* geometry_,double** pumps_){
    geometry=geometry_;
    pumps=pumps_;
    indice_x_Richards=new bool[geometry->nX];
    richards_evolutive_time=new RichardsEvolutiveTime[geometry->nX];
    for(size_t ix=0;ix<geometry->nX;++ix){
      richards_evolutive_time[ix].init(ix,geometry,pumps[ix]);
    }
  }

  void
  AllVerticalRichards::update_indice_x_Richards(){
    for(size_t ix=0;ix<geometry->nX;++ix){
      indice_x_Richards[ix]=(error_x[ix]>max_error_x);
    }
  }


  void
  AllVerticalRichards::init_indice_x_Richards(){
    for(size_t ix=0;ix<geometry->nX;++ix){
      indice_x_Richards[ix]=true;
    }
  }


  void
  AllVerticalRichards::run(){
    Debug::debug_Thomas=false;

    size_t i_left,i_middle,i_right;

    #ifdef MYDEBUG
    cout<<"   [AllVerticalRichards::run] start"<<endl;
    #endif

    for(size_t ix=0;ix<geometry->nX;++ix){

      P[ix][geometry->nZ[ix]-1]=Psoil[ix];
      if(indice_x_Richards[ix]){
        RichardsEvolutiveTime& r=richards_evolutive_time[ix];
        r.dt=dt;
        r.init_P=init_P[ix];
        r.previous_P=previous_P[ix];
        r.P=P[ix];


        //Update Psoil
  /*      r.P[geometry->nZ[ix]-1]=Psoil[ix];
        r.init_P[geometry->nZ[ix]-1]=Psoil[ix];
        r.previous_P[geometry->nZ[ix]-1]=Psoil[ix];*/



        if(ix==0){
          i_left=0;
          i_middle=1;
          i_right=2;
        }
        else if(ix==geometry->nX-1){
          i_left=geometry->nX-3;
          i_middle=geometry->nX-2;
          i_right=geometry->nX-1;
        }
        else{
          i_left=ix-1;
          i_middle=ix;
          i_right=ix+1;
        }
        r.hydr_left=hydr[i_left];
        r.l_left=l[i_left];
        r.Pl_left=Pl[i_left];
        r.hydr_middle=hydr[i_middle];
        r.l_middle=l[i_middle];
        r.Pl_middle=Pl[i_middle];
        r.hydr_right=hydr[i_right];
        r.l_right=l[i_right];
        r.Pl_right=Pl[i_right];
        r.run();
        P[ix]=r.P;
      }
    }
    //[If Paraller wait for sync]
    has_converged=true;
    for(size_t ix=0;ix<geometry->nX;++ix){
      if(!richards_evolutive_time[ix].has_converged){
        has_converged=false;
      }
    }
    if(has_converged) compute_hsat();

    #ifdef MYDEBUG
    cout<<"   [AllVerticalRichards::run] ";
    if(not has_converged) cout<<"not";
    cout<<"converge"<<endl<<"   [AllVerticalRichards::run] end"<<endl;
    #endif

  }

  void
  AllVerticalRichards::compute_hsat(){
    for(size_t ix=0;ix<geometry->nX;++ix){
      double* Px=previous_P[ix];
      size_t iz=0;
      for(;iz<geometry->nZ[ix] and Px[iz]>Psat;++iz);
      if(iz>0 and Px[iz]<=Psat){
        hsat[ix]=(Psat-Px[iz-1])*geometry->dZ[ix]/(Px[iz]-Px[iz-1])+geometry->Z[ix][iz-1];
      }
      else{
        hsat[ix]=(iz==0)?geometry->hbot[ix]:geometry->hsoil[ix];
      }
    }
  }
}