123456789101112131415161718192021222324252627282930313233343536373839404142434445464748 |
- #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){
- //TODO :: add rain
- hov[ix]=init_hov[ix]-dt*flux_soil[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;
- #ifdef MYDEBUG
- if(Debug::ix==ix){
- cout<<" [Overland::compute_flux_soil] p1 = "<<p1<<endl;
- cout<<" [Overland::compute_flux_soil] p2 = "<<p2<<endl;
- cout<<" [Overland::compute_flux_soil] flux = "<<flux<<endl;
- }
- #endif
- }
- }
- }
|