module soil_vegetation_module

contains

!===================================================================================!

subroutine tau_soil_veg( s, par , tau_nc , tau_c) ! , fract )
    !
    use params
    use spaceparams
    use vegetation_module
    !
    implicit none
    !
    type(parameters)                      :: par
    type(spacepars),target                :: s
    !type(veggie), dimension(:), pointer  :: veg
    real*8, dimension(par%nx+1,par%ny+1), intent(inout) :: tau_nc
    real*8, dimension(par%nx+1,par%ny+1), intent(inout) :: tau_c
    real*8, dimension(par%nx+1,par%ny+1) :: psi2C     ! [] constant
    real*8, dimension(par%nx+1,par%ny+1) :: psi2NC    ! [] constant
    real*8                              :: da        ! [m] average size of particel aggregates
    real*8                              :: cu        ! [N/m2] soil undrained shear strenght 
    real*8                              :: tr        ! [N/m2] root tensile strength
    real*8                              :: RAR       ! [-] root surface area ratio
    real*8                              :: mrho      ! [kg/m3] soil mean density
    real*8                              :: delta     ! [-] submerged weight
    real*8                              :: Cr        ! [N/m2] "root cohesion"
    real*8                              :: Cf        ! [N/m2] soil cohesion
    integer                             :: i,j,jg
    !
    !
    ! Soil characteristics
    !
    cu   = 20.0E03
    !mrho = bc%settings%rhofrac(1)*fract(1)+bc%settings%rhofrac(2)*fract(2)
    mrho = par%rhos
    da   = 0.01d0
    !
    ! Vegetation characteristics
    ! TODO: all this value have to come from some external files (check file for vegetation)
    RAR = 0.01d0
    tr  = 1500.0E03
    !
    ! Determination fo the constant psi2
    !
    delta  = (mrho-par%rho)/par%rho
    Cf     = 0.035d0*cu
    psi2C  = tau_c/(delta*par%g*da+1.0d0/par%rho*0.6d0*Cf)
    psi2NC = tau_nc/(delta*par%g*da+1.0d0/par%rho*0.6d0*Cf)
    !
    ! Determination of "root cohesion"
    !
    Cr = 1.2*RAR*tr
    !
    ! Determination of modified critical shear stress
    !
    do j=1,par%ny+1
        do i=1,par%nx+1
            if (par%soilveg==1 .and. s%vegtype(i,j)==1) then
               tau_c(i,j)  = psi2C(i,j)*(delta*par%g*da+1/par%rho*(0.6d0*Cf+0.6d0*Cr))
               tau_nc(i,j) = psi2NC(i,j)*(delta*par%g*da+1/par%rho*(0.6d0*Cf+0.6d0*Cr))
            end if
        end do
    end do   
    !
end subroutine tau_soil_veg

!===================================================================================!

end module soil_vegetation_module