RichardsFastSlow/src/kernel/overland.cpp

#include "overland.hpp"
#include "debug.hpp"

namespace Kernel{

  void Overland::run(){
    #ifdef MYDEBUG
    cout<<"   [Overland::run] start"<<endl;
    #endif
    compute_flux_soil();
    total_error=0;
    for(size_t ix=0;ix<geometry->nX;++ix){
      hov[ix]=init_hov[ix]-dt*flux_soil[ix];
      //TODO :: add rain
      if(ix>89){
        hov[ix]+=dt*0.0001;
      }
    }
    apply_flow();
    for(size_t ix=0;ix<geometry->nX;++ix){
      double e=abs(hov[ix]-previous_hov[ix])/n1_init_hov;
      error_x[ix]+=e;
      total_error+=e;

      Psoil[ix]=Physics::Psoil(geometry->hsoil[ix],hov[ix]);
    }
    #ifdef MYDEBUG
    cout<<"   [Overland::run] error : "<<total_error<<endl;
    #endif
    #ifdef MYDEBUG
    cout<<"   [Overland::run] end"<<endl;
    #endif
  }

  void Overland::compute_flux_soil(){
    for(size_t ix=0;ix<geometry->nX;++ix){
      double p1=P[ix][geometry->nZ[ix]-1];
      double p2=P[ix][geometry->nZ[ix]-2];

      double flux=Physics::k0/2*(Physics::kr(p1)+Physics::kr(p2))*((p1-p2)/(geometry->dZ[ix]*Physics::rho*Physics::g)+1);
      if(ix==0) flux-=Physics::k0*Physics::kr(Pl[0]+Physics::rho*Physics::g*(l[0]-geometry->hsoil[0]))*(hydr[1]-hydr[0])*geometry->dhsoil[0]/geometry->dX;
      else if(ix==geometry->nX-1) flux-=Physics::k0*Physics::kr(Pl[ix]+Physics::rho*Physics::g*(l[ix]-geometry->hsoil[ix]))*(hydr[ix]-hydr[ix-1])*geometry->dhsoil[ix]/geometry->dX;
      else flux-=Physics::k0*Physics::kr(Pl[ix]+Physics::rho*Physics::g*(l[ix]-geometry->hsoil[ix]))*(hydr[ix+1]-hydr[ix-1])/2*geometry->dhsoil[ix]/geometry->dX;
      flux_soil[ix]=flux;

    }
  }

  void Overland::apply_flow(){
    bool overflow=true;
    while(overflow){
      overflow=false;
      runoff_to_local_min();
      basin_transfer();
    }
  }

  void Overland::runoff_to_local_min(){
  //  cout<<"      [Overland::runoff_to_local_min]"<<endl;
    for(size_t ix=0;ix<geometry->nX;++ix){
      double dh=hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff;
      if(dh>0){
        double dh_left,dh_right;
        switch(geometry->runoff_directions[ix]){
          case Left:
          dh_left=dh;
          dh_right=0;
          break;
          case Right:
          dh_right=dh;
          dh_left=0;
          break;
          case Both:
          dh_left=dh/2;
          dh_right=dh/2;
          break;
          default:
          dh_left=0;
          dh_right=0;
          break;
        }
        size_t jx=ix;
        while(dh_left>0){
          --jx;
          if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Right){
            hov[jx]+=dh_left;
            hov[ix]-=dh_left;
            break;
          }
          else{
            double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_left);
            hov[jx]+=temp;
            hov[ix]-=temp;
            dh_left-=temp;
          }
        }
        jx=ix;
        while(dh_right>0){
          ++jx;
          if(geometry->runoff_directions[jx]==None or geometry->runoff_directions[jx]==Left){
            hov[jx]+=dh_right;
            hov[ix]-=dh_right;
            break;
          }
          else{
            double temp=min(max(Physics::minimal_height_to_runoff-hov[jx]+geometry->hsoil[jx],0.0),dh_right);
            hov[jx]+=temp;
            hov[ix]-=temp;
            dh_right-=temp;
          }
        }

      }
    }
  /*cout<<"------------------- after runoff ----------------"<<endl;
    for(size_t ix=0;ix<100;++ix){
      cout<<ix<<" : "<<max(0.0,hov[ix]-geometry->hsoil[ix]-Physics::minimal_height_to_runoff)<<endl;;
    }*/
  }

  void Overland::basin_transfer(){
    double vol_tot=0;
    for(size_t ix=0;ix<geometry->nX;++ix){
      vol_tot+=max(hov[ix]-geometry->hsoil[ix],0.0);
    }
    //TODO :: Remove guard
    int nmax=100;
    do{
      --nmax;
      geometry->root_basin->compute_vflow(geometry->hsoil,hov);
      geometry->root_basin->repartition(geometry->hsoil,hov,geometry->runoff_directions);
    }while(geometry->root_basin->overflow and nmax>0);
    if(nmax==0){
      cerr<<"Too many loop"<<endl;
      Debug::pause();//exit(-1);
    }
    geometry->root_basin->compute_vflow(geometry->hsoil,hov);
    double vol_tot_mid=0;
    for(size_t ix=0;ix<geometry->nX;++ix){
      vol_tot_mid+=max(hov[ix]-geometry->hsoil[ix],0.0);
    }
    geometry->root_basin->flatten(geometry->hsoil,hov);

    double vol_tot_new=0;
    for(size_t ix=0;ix<geometry->nX;++ix){
      vol_tot_new+=max(hov[ix]-geometry->hsoil[ix],0.0);
    }
  }
}