!! Set soil properties 
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  1.3 - 3rd October 2016  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 06/Dec/2013 | Original code |
! | 1.1      | 25/Feb/2016 | added parameters for van Genuchten WRC |
! | 1.2      | 16/Jun/2016 | added parameters for Curve Number |
! | 1.3      | 16/Jun/2016 | added function DepthToDiv |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! Define structure of soil column and soil retention curves
!
MODULE SoilProperties        
			

! Modules used: 

USE DataTypeSizes, ONLY : &
! Imported Type Definitions:
short, float, double   

USE LogLib, ONLY: &
!Imported routines:
Catch

USE StringManipulation, ONLY: &
!Imported routines:
ToString

IMPLICIT NONE 
! Global (i.e. public) Declarations:

TYPE SoilType
  !hydrological properties
  REAL (KIND = double) :: ksat !! saturated hydraulic conductivity [m/s]
  REAL (KIND = double) :: thetas !! saturated volumetric water content [m3/m3]
  REAL (KIND = double) :: thetar !! residual volumetric water content [m3/m3]
  REAL (KIND = double) :: wp !! wilting point [m3/m3]
  REAL (KIND = double) :: fc !! field capacity [m3/m3]
  REAL (KIND = double) :: psic !! air entry value [m]
  REAL (KIND = double) :: pp !! tortuosity index, pre interaction parameter [-]
  !required by Green-Ampt
  REAL (KIND = double) :: phy !!suction head across the wetting front [m]
  REAL (KIND = double) :: smax !!maximum soil surface storage [m]
  !required by Richards with Brooks and Corey WRC
  REAL (KIND = double) :: psdi !! Brooks and Corey pore size distribution index [-]
  !required by Richards with van Genuchten WRC
  REAL (KIND = double) :: m
  REAL (KIND = double) :: n
  REAL (KIND = double) :: kx !!saturated K of soil "matrix" [m/s]
  !Curve Number parameters
  REAL (KIND = double) :: c !!Curve Number initial abstraction ratio
  REAL (KIND = double) :: s0 !!Storativity, default = 254 [mm]
  REAL (KIND = double) :: cn !!Curve Number value
  
END TYPE SoilType

TYPE SoilLayer
  !geometry
  REAL (KIND = double) :: thickness !! layer thickness [m]
  !soil type
  INTEGER (KIND = short) :: stype !!soil type
  !state variable
  REAL (KIND = double) :: theta !! actual volumetric water content [m3/m3]
  REAL (KIND = double) :: kact !!actual hydraulic conductivity [m/s]
  REAL (KIND = double) :: psi !!actual matric potential [m] 
END TYPE SoilLayer

!! convention: first layer is close to ground surface
TYPE SoilColumn
  INTEGER :: types !! number of soil types
  INTEGER :: divs !! number of subdivisions in soil column
  TYPE(SoilLayer), ALLOCATABLE :: div (:)
  TYPE(SoilType),  ALLOCATABLE :: typ (:)
END TYPE SoilColumn 

!public routines:
PUBLIC :: UnsHydCond
PUBLIC :: Psi
PUBLIC :: DepthToDiv
                
!=======
CONTAINS
!=======
! Define procedures contained in this module. 

!==============================================================================
!| Description:
! Compute hydraulic conductivity of partially saturated soil (m/s)
FUNCTION UnsHydCond &
!
 (ksat, theta, thetas, thetar, psdi) &
 !
 RESULT (k)

IMPLICIT NONE

!Arguments with intent in
REAL (KIND = float), INTENT(IN) :: ksat !!saturated hydraulic conductivity (m/s)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: psdi !!Brooks & Corey pore size distribution index (-)


!local declarations:
REAL (KIND = float) :: k


!-------------------------end of declarations----------------------------------

k = ksat * ( (theta - thetar) / (thetas - thetar) ) ** (2./psdi + 3.)


RETURN

END FUNCTION UnsHydCond


!==============================================================================
!| Description:
! Compute matric potential (m)
FUNCTION Psi &
!
 (psic, theta, thetas, thetar, psdi) &
!
 RESULT (h)

IMPLICIT NONE

!Arguments with intent in
REAL (KIND = float), INTENT(IN) :: psic !!bubbling pressure (m)
REAL (KIND = float), INTENT(IN) :: theta !!volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetas !!saturated volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: thetar !!residual volumetric water content (m3/m3)
REAL (KIND = float), INTENT(IN) :: psdi !!Brooks & Corey pore size distribution index (-)


!local declarations:
REAL (KIND = float) :: h

!-------------------------end of declarations----------------------------------


 h = psic / ( (theta - thetar) / (thetas - thetar) ) ** (1./psdi)

RETURN

 END FUNCTION Psi



!==============================================================================
!| Description:
! Return which division corresponds to given depth
 FUNCTION DepthToDiv &
!
 (column, depth) &
!
 RESULT (d)

IMPLICIT NONE

!Arguments with intent in
TYPE(SoilColumn), INTENT(IN)  :: column
REAL (KIND = float), INTENT(IN) :: depth

!local declarations:
INTEGER :: d
LOGICAL :: founddiv
REAL :: depth_computed

!-------------------------end of declarations----------------------------------
founddiv = .FALSE.
depth_computed = 0.
DO d = 1, column % divs
  depth_computed = depth_computed + column % div (d) % thickness
  IF (depth_computed == depth) THEN
    founddiv = .TRUE.
    RETURN
  END IF
END DO

IF (.NOT. founddiv) THEN
  CALL Catch ('error', 'SoilProperties', &
       'division not found corresponding to depth: ' , &
        argument = ToString(depth))
     
  
END IF

END FUNCTION DepthToDiv

END MODULE SoilProperties