!> algebraic functions for the porous media flow model by M. Kuraz
!> data parameters
module porous_fnc
  use spline_rec
  
  implicit none

  type, public :: soilpar_str
     real :: alpha  ! Van Genuchten
     real :: n      ! Van Genuchten
     real :: m      ! Van Genuchten
     real :: Thr    ! theta_r
     real :: Ths    ! theta_s
     real :: Ss     ! storativity
     real :: Ks     ! saturated conductivity

     real :: capacity_max ! maximal capacity for cubic interpolation
     real :: h_max        ! position of the achieved  maximal capacity
     real :: h_max2       ! limit for theconductivity modification

     real :: capa_a    ! cubic interpolation of the conductivity, coeff = a 
     real :: capa_b    ! cubic interpolation of the conductivity, coeff = b 
     real :: capa_d    ! cubic interpolation of the conductivity, coeff = d 

     real :: cond_a    ! cubic interpolation of the conductivity, coeff = a 
     real :: cond_b    ! cubic interpolation of the conductivity, coeff = b 
     real :: cond_d    ! cubic interpolation of the conductivity, coeff = d 
     real :: cond_e    ! cubic interpolation of the conductivity, coeff = e .. correction

     real :: num_vol
     real :: h_int
     
     real, dimension(:), pointer :: volume
     real, dimension(:,:), pointer :: rec_spline
  end type soilpar_str
  
  logical, public :: USE_TRACY=.false.

  !type(soilpar_str), private, dimension(:), allocatable :: soilpar
  type(soilpar_str), dimension(:), allocatable :: soilpar

  !> number of materials
  integer :: number_soilpar 

  !> dynamic viscosity
  real, parameter, private :: mu=1.62e-7
  !> density of water
  real, parameter, private :: rho=1000.0
  !> "configuration of the void space" exact value definition should be improved
  !real, parameter, private :: beta=1.0
  ! VD
  real, parameter, private :: beta=0.0

  !> cubic interpolation switch on/off
  logical :: cubic_interpol
  
  !> spline interpolation of the water contents
  logical :: water_cubic_interpol
  
  public :: init_porous_coeffs
  !> van Genucheten soil hydraulic functions
  public :: vangen, conduct, capacity
  !> Forcheimer update for the hydraulic conductivity
  public :: forch_conduct
  !> Tracy's model for the soil hydraulic functions
  public :: gardner_cond, gardner_wc, gardner_cap

contains

  !> assigns soil parameters and allocates structures,
  !> at this moment this function is just simple initialization of values for the case study of dam seepage, can be always improved:)
  subroutine init_porous_coeffs()
    implicit none
    integer :: layers

    ! "correction" of the conductivity and capacity
    cubic_interpol = .false.  

    ! spline interpolation of the water contents
    water_cubic_interpol = .false.
    
    layers = 4 ! in current case study we have 3 different materials

    number_soilpar = layers
    
    allocate(soilpar(layers))

!!!!!!!!!!!!!!!!!!!!!!
    !units: m, days
!!!!!!!!!!!!!!!!!!!!!

    !layer 1 = gravel
    ! original values
    !soilpar(1)%alpha = 100.0
    !soilpar(1)%n = 2.0
    !soilpar(1)%m = 0.626

    ! New values by Michal
    if (USE_TRACY) then
      soilpar(1)%alpha = 0.1
    else
      soilpar(1)%alpha = 2.0
    end if

    soilpar(1)%n = 1.41
    soilpar(1)%m = 0.291


    ! New values by Michal (2)
    !soilpar(1)%alpha = 0.8  
    !soilpar(1)%n = 1.2      
    !soilpar(1)%m = 0.167    

    if (USE_TRACY) then
      soilpar(1)%Ks = 1.0
      soilpar(1)%ths = 0.5
      soilpar(1)%thr = 0.0
    else
      soilpar(1)%Ks = 7.128
      soilpar(1)%ths = 0.43
      soilpar(1)%thr = 0.01   ! 0.
    end if
    
    !soilpar(1)%ths = 0.43  !0.5
    !soilpar(1)%thr = 0.01   ! 0.
    soilpar(1)%Ss = 1.0E-02

    
    !layer 2 = clay
    soilpar(2)%alpha = 0.8  !0.8   !!0.8
    soilpar(2)%n = 1.05  !1.2    !n \in (1.05, 2.5)  !! orig = 1.2
    soilpar(2)%m = 0.167
    soilpar(2)%Ks = 0.048
    soilpar(2)%ths = 0.38
    soilpar(2)%thr = 0.06
    soilpar(2)%Ss = 1.0E-02

    !layer 3 = silt clay
    soilpar(3)%alpha = 2.0
    soilpar(3)%n = 1.41
    soilpar(3)%m = 0.291
    soilpar(3)%Ks = 0.108
    soilpar(3)%ths = 0.45
    soilpar(3)%thr = 0.067
    soilpar(3)%Ss = 1.0E-02

    !layer 4 = ../porous/valcovazk.ini
    soilpar(4)%alpha = 0.8
    soilpar(4)%n = 1.2
    soilpar(4)%m = 0.167
    soilpar(4)%Ks = 2.77777777784352E-05
    soilpar(4)%ths = 0.55
    soilpar(4)%thr = 0.0
    soilpar(4)%Ss = 1.0E-03

    !soilpar(4)%alpha = 0.8  !0.8   !!0.8
    !soilpar(4)%n = 1.05  !1.2    !n \in (1.05, 2.5)  !! orig = 1.2
    !soilpar(4)%m = 0.167
    soilpar(4)%Ks = 0.048
    !soilpar(4)%ths = 0.38
    !soilpar(4)%thr = 0.06
    !soilpar(4)%Ss = 1.0E-02

    
    ! ! REALISTIC values
    ! ! coarse sand
    ! soilpar(1)%alpha = 200.
    ! soilpar(1)%n = 2.6 
    ! soilpar(1)%m = 0.615
    ! soilpar(1)%Ks = 1.1574074074074E-06
    ! soilpar(1)%ths = 0.55
    ! soilpar(1)%thr = 0.0
    ! soilpar(1)%Ss = 1.0E-02


    ! ! clay
    ! soilpar(2)%alpha = 2.
    ! soilpar(2)%n = 1.2 
    ! soilpar(2)%m = 0.167
    ! soilpar(2)%Ks = 1.5E-09
    ! soilpar(2)%ths = 0.45
    ! soilpar(2)%thr = 0.05
    ! soilpar(2)%Ss = 1.0E-03


    ! ! silt clay
    ! soilpar(3)%alpha = 1.
    ! soilpar(3)%n = 1.6 
    ! soilpar(3)%m = 0.375
    ! soilpar(3)%Ks = 1.0E-07
    ! soilpar(3)%ths = 0.48
    ! soilpar(3)%thr = 0.01
    ! soilpar(3)%Ss = 5.0E-03


    
    !print*,'NO storativity  !!!!!!!!!!!!!!!!!!'
    !soilpar(1:3)%Ss = 0. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1

    ! New values by Michal
    ! soilpar(1:3)%alpha = 2.0
    ! soilpar(1:3)%n = 1.41
    ! soilpar(1:3)%m = 0.291
    ! soilpar(1:3)%Ks = 7.128
    ! soilpar(1:3)%ths = 0.43
    ! soilpar(1:3)%thr = 0.0
    ! soilpar(1:3)%Ss = 1.0E-02

    !layer 2 = clay
    !soilpar(1:3)%alpha = 0.8
    !soilpar(1:3)%n = 1.2
    !soilpar(1:3)%m = 0.167
    !soilpar(1:3)%Ks = 0.048
    !soilpar(1:3)%ths = 0.38
    !soilpar(1:3)%thr = 0.06
    !soilpar(1:3)%Ss = 1.0E-02

    !soilpar(1:3)%alpha = 2.0
    ! soilpar(1:3)%n = 1.41
    ! soilpar(1:3)%m = 0.291
    ! soilpar(1:3)%Ks = 0.108
    ! soilpar(1:3)%ths = 0.45
    ! soilpar(1:3)%thr = 0.067
    ! soilpar(1:3)%Ss = 1.0E-02

    ! soilpar(1:3)%alpha = 2.0
    ! soilpar(1:3)%n = 4.
    ! soilpar(1:3)%m = 2.
    ! soilpar(1:3)%Ks = 1.0
    ! soilpar(1:3)%ths = 0.5
    ! soilpar(1:3)%thr = 0.1
    ! soilpar(1:3)%Ss = 1.0E-00

  end subroutine init_porous_coeffs


  !> \brief Van Genuchten relation \f$ \theta = f(pressure) \f$
  !>  \f$ \theta_e = \frac{1}{(1+(\alpha*h)^n)^m} \f$
  !> water content is considered as absolute value not the relative one \n
  !> see \f$ \theta_e = \frac{\theta - \theta_r}{\theta_s-\theta_r} \f$
  !>
  function vangen(u, layer, z) result(theta)
    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting water content
    real :: theta

    real :: a,n,m, theta_e, h

    h = u - z


    a = soilpar(layer)%alpha
    n = soilpar(layer)%n
    m = soilpar(layer)%m


    if (h >=0.0) then
       theta = soilpar(layer)%Ths
       RETURN
    else
       theta_e = 1/(1+(a*(abs(h)))**n)**m
       theta = theta_e*( soilpar(layer)%Ths - soilpar(layer)%Thr ) +soilpar(layer)%Thr
    end if

  end function vangen

  !> \brief so-called retention water capacity, it is a derivative of the retention curve function
  !> \f$ E(h) = C(h) + \frac{\theta(h)}{\theta_s}S_s \f$
  !> where
  !> \f$ C(h) = \left\{ \begin{array}{l l}\frac{m n \alpha  (-h \alpha )^{-1+n}}{\left(1+(-h \alpha )^n\right)^{1+m}}(\theta_s - \theta_r),
  !> & \quad \mbox{$\forall$ $h \in (-\infty, 0 )$}\\ 0, & \quad \mbox{$\forall$ $h \in \langle 0, + \infty )$}\\ \end{array} \right. \f$
  !> and
  !> \f$ \theta(h) = \left\{ \begin{array}{l l} \frac{\theta_s -\theta_r}{(1+(-\alpha h)^n_{vg})^m_{vg}} + \theta_r,
  !> & \quad \mbox{$\forall$ $h \in (-\infty, 0 )$}\\ \theta_S, & \quad \mbox{$\forall$ $h \in \langle 0, + \infty )$}\\ \end{array} \right. \f$
  !>
  function capacity(u, layer, z) result(E)

    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting system capacity (elasticity)
    real :: E, r

    real :: C, a, m, n, tr, ts, h

    h = u - z

    if (h < 0) then
       a = soilpar(layer)%alpha
       n = soilpar(layer)%n
       m = soilpar(layer)%m
       tr = soilpar(layer)%Thr
       ts = soilpar(layer)%Ths
       C = a*m*n*(-tr + ts)*(-(a*h))**(-1 + n)*(1 + (-(a*h))**n)**(-1 - m)
    else
       E = soilpar(layer)%Ss
       RETURN
    end if

    E = C + vangen(u, layer, z)/soilpar(layer)%Ths*soilpar(layer)%Ss

    ! modification by the cubic interplation in vicinity  h < 0. 
    if( cubic_interpol) then
       if( h< 0. .and. h >=  soilpar(layer)%h_max) then
          r = h /  soilpar(layer)%h_max
          
          E = ( soilpar(layer)%Ss - soilpar(layer)%capacity_max )  &
               * r*r * (2 *r - 3) + soilpar(layer)%Ss

          E = r*r *( soilpar(layer)%capa_a * r +  soilpar(layer)%capa_b) +  soilpar(layer)%capa_d 
               !+ soilpar(layer)%cond_e * r * r* (1-r)*(1-r)
       endif
          
    end if
    ! VD
    !E = soilpar(layer)%Ss
    !write(*,'(a8, 20es12.4)') 'params:', soilpar(layer)%Ss, E , E/ soilpar(layer)%Ss
    !write(*,'(a8, 20es12.4)') 'params:', soilpar(layer)%Ss, a, n, m, tr, ts, soilpar(layer)%Ks

  end function capacity


  !> \brief Mualem's function for unsaturated hydraulic conductivity with van Genuchten's water content substitution
  !> \f$   K(h) = \left\{ \begin{array}{l l} K_s\frac{\left( 1- (-\alpha h)^{n_{vg}m_{vg}} \left( 1+ (-\alpha h)^{n_{vg}}
  !> \right)^{-m_{vg}} \right)^2}{\left(1+(-\alpha h)^{n_{vg}} \right)^{\frac{m_{vg}}{2}}},  &  \mbox{$\forall$
  !> $h \in$ $(-\infty,0)$}\\ K_s, & \mbox{$\forall$   $h \in$ $\langle 0, +\infty)$}\\ \end{array} \right. \f$
  !>
  function conduct(u, layer, z) result(K)

    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting hydraulic conductivity
    real :: K
    real :: a,n,m,h
    real :: K1, ah, r

    h = u - z

    if (h >  0) then
       K = soilpar(layer)%Ks
    else
       a = soilpar(layer)%alpha
       n = soilpar(layer)%n
       m = soilpar(layer)%m

       K =  (1 - (-(a*h))**(m*n)/(1 + (-(a*h))**n)**m)**2/(1 + (-(a*h))**n)**(m/2.0) * soilpar(layer)%Ks
       !ah = -a*h

       !print*,'#DE#:',( 1 + (ah)**n)**(m)
       !print*,'#DE#:',( 1 + (ah)**n )**(m/2), ah, n*m
       !print*,'????:',  ( 1  + ah**(n*m) / ( 1 + (ah)**n)**(m) )**2

       !K1 = ( 1 - ah**(n*m) / ( 1 + ah**n)**m )**2 / ( 1 + ah**n )**(m/2) * soilpar(layer)%Ks

       !if( K /soilpar(layer)%Ks  < 1E-10) &
       !     write(*,'(a8, 26es12.4)') ' K = ', K, K1, soilpar(layer)%Ks,  K/  soilpar(layer)%Ks, &
       !     h, u, z
    end if

    ! modification by the cubic interplation in vicinity  h < 0. 
    if( cubic_interpol) then
       if( h< 0. .and. h >=  soilpar(layer)%h_max2) then
         r = h /  soilpar(layer)%h_max2

         K = r*r *( soilpar(layer)%cond_a * r +  soilpar(layer)%cond_b) +  soilpar(layer)%cond_d &
              + soilpar(layer)%cond_e * r * r* (1-r)*(1-r)

         !write(*,'(a8, 30es12.4)') &
         !     'cdeek:',soilpar(layer)%cond_e, r, soilpar(layer)%cond_e * r * r* (1-r)*(1-r), K, &
         !     soilpar(layer)%cond_e * r * r* (1-r)*(1-r) / K, soilpar(layer)%h_max2

      end if
    end if

    
    ! VD
    !K = soilpar(layer)%Ks

  end function conduct

  function forch_conduct(u, gradu, layer, z) result(kappa)
    implicit none
    !> solution
    real, intent(in) :: u
    !> L2 norm of the solution gradient
    real, intent(in) :: gradu
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z

    !> resulting kappa conductivity for Forchheimer equation
    real :: kappa

    real :: K

    K = conduct(u, layer, z)


    ! ERROR ?, exchanged a_0 and a_1 !!!
    !kappa = 2.0/(1.0/K + sqrt( (rho*beta/(mu*K))*(rho*beta/(mu*K)) + 4.0/K*gradu))

    ! better computer arithmetic, works for K= 0
    !kappa = 2.0 * K /(1.0 + sqrt( (rho*beta/(mu))*(rho*beta/(mu)) + 4.0*K*gradu))

    ! VD CORRECTED
    kappa = 2.0 * K /(1.0 + sqrt(1 +  4. * (rho*beta/ mu)*K *gradu))

  end function forch_conduct


  !> spline reconstruction of the numerical integration of \f$ \theta(h) \f$ for
  !> the water content
  subroutine porous_water_content_spline( ) 
    implicit none
    real, dimension(:,:), allocatable :: hi
    integer :: layer, j, nn, n8
    real :: h, h_max, val_max, h_left, h_right, h_infty, h_step
    real :: hL, hR
    real :: val, sum

    
    
    do layer=1, number_soilpar

       ! number of precomputed values to be interpolated
       nn = soilpar(layer)%num_vol + 1

       ! values of the spline reconstruction
       allocate(soilpar(layer)%rec_spline(1:3, 1: nn+1) )

       ! x- values
       allocate(hi(1:2, 0: nn+2) )
       do j=1, nn
          hi(1, j) = (j-1.0) / soilpar(layer)%h_int 
          hi(2, j) = soilpar(layer)%volume(j-1)
       enddo
       
       call spline_reconstruct(nn, hi(1, 0:nn+1), hi(2, 0:nn+1), &
            soilpar(layer)%rec_spline(1, 1: nn+1), &
            soilpar(layer)%rec_spline(2, 1: nn+1), &
            soilpar(layer)%rec_spline(3, 1: nn+1), layer )

       !do j=1, nn
       !   write(30+layer, *) hi(1:2, j+1), soilpar(layer)%volume(j-1) 
       !enddo

       !stop "water integ spline: e3t33y"
       
       deallocate( hi) 
    enddo

    water_cubic_interpol = .true.

    !stop "water integ spline:"
    
  end subroutine porous_water_content_spline

  
  !> numerical integration of \f$ \theta(h) \f$ for the water content
  subroutine porous_water_content_integ(icase ) 
    implicit none
    integer, intent(in) :: icase
    integer :: layer, j, nn, n8
    real :: h, h_max, val_max, h_left, h_right, h_infty, h_step
    real :: hL, hR
    real :: val, sum

    !h_infty = -1E+5  ! limit value approxmating infty
    h_left = -100.   ! starting value, sufficiently far from the real region
    h_right = 0.

    if(icase == 5) then
       h_left = -50.
    endif

    !print*, 'H_left = ', h_left
       
    
    !nn = 100
    nn = 500
    !nn = 1500
    !nn = 2500
    !nn = 5000
    !n8 = 100
    
    do layer=1, number_soilpar
       h_step = (h_right - h_left) / nn

       soilpar(layer)%h_int = nn/  h_left 

       soilpar(layer)%num_vol = nn 
       allocate(soilpar(layer)%volume(0: soilpar(layer)%num_vol) )

       ! computing the water contents from h_infty to h_left
       hL = h_left
       hR = h_left + h_step

       sum = 0.  ! initial approximation

       soilpar(layer)%volume(nn) = sum

       do j = 1, nn
          call integ_volume_content(layer, hL, hR, val )
          sum = sum + val

          soilpar(layer)%volume(nn-j) = sum

          hL = hR
          hR = hL + h_step

       enddo

       !write(*,'(a8, i5, 3es12.4, i5)') 'hd383:', layer, val, sum, h_step, j
       !do j=0, nn
       !   hL = h_left + 1. * j * h_step
       !   write(30 + layer, *) hL,  soilpar(layer)%volume(j)
       !enddo
    enddo
    !stop "water integ:"
    
  end subroutine porous_water_content_integ


  

  !> integration of the function \f$ \theta(h) \f$ over \f$ (h_L, h_R) \f$
  subroutine integ_volume_content(layer, hL, hR, val )
    integer, intent(in) :: layer   ! material id
    real, intent(in) :: hL, hR     ! size of the interval
    real, intent(inout) :: val     ! value of the integral
    !!!type(Gauss_rule), pointer :: G_rule
    integer :: Gnum, i
    real :: h, theta
    real :: weights(1:3), lambda(1:3)
    
    Gnum = 3
    weights(1) = 5./18
    weights(2) = 8./18
    weights(3) = 5./18
    lambda(1) = (1.-0.77459666924148337704)/2
    lambda(2) = 0.5
    lambda(3) = 1.77459666924148337704/2

    
    !!G_rule => state%space%G_rule(Gnum)

    val = 0.
    do i=1,Gnum
       h = hL + lambda(i) * (hR - hL)
       theta = capacity(h, layer, 0.)

       val = val + weights(i) * theta
       !write(20 + layer, *) h, theta
       !write(22, *) i, h, theta, val * (hR - hL)
    enddo
    
    val = val *(hR - hL)

    !write(23, *) hL, val
    !write(23, *) hR, val
    !stop "738733"
    
  end subroutine integ_volume_content
  
  !> evaluation of maximal values for the cubic interpolation near h < 0
  subroutine porous_coeffs_cubic_interpol(param1)
    implicit none
    integer :: i, j, nn
    real :: h, h_max, val, val_max, h_left, h_right,param1
    real :: kappa_der, kappa_M, kappa_0, hh, A, B, alpha

    
    nn = 1000
    do i=1, number_soilpar
       h_left = -1.  ! initial guess of the interval containing max_value
       h_right = 0. 
       
       
       val_max = 0.
       do j=0, nn
          h = h_right - 1.*j/nn * (h_right - h_left)
          val = capacity(h, i, 0.)
          if(val > val_max) then
             val_max = val
             h_max = h
          endif
       enddo
       if( val == val_max) then
          print*,' maximum at boundary node, probably did not find (due to monotonicity)'
       !   stop
       endif

       ! capacity
       soilpar(i)%capacity_max = val_max
       soilpar(i)%h_max = h_max

       ! new variants with variable h_max
       !!h_max = -0.5
       !h_max = -0.05
       !!h_max = -0.025
       !h_max = -0.1
       h_max = -0.01
       !h_max = -0.001
       !h_max = -0.0001
       h_max = -param1 
       print*,'h-Max = ', h_max, soilpar(i)%h_max
       
       soilpar(i)%h_max = h_max

       ! conductivity in the left and right
       kappa_0 = capacity(h_right, i, 0.)
       kappa_M = capacity(h_max, i, 0.)

       ! derivative of the condutivity in the left
       hh = abs( (h_right - h_left) / 200.) 
       kappa_der = (capacity(h_max+hh, i, 0.) - capacity(h_max-hh, i, 0.) ) / (2*hh)

       print*,i,'kappa_der =',kappa_der, h_max

       soilpar(i)%capa_a = -2.*(kappa_M - kappa_0) + kappa_der * h_max
       soilpar(i)%capa_b =  3.*(kappa_M - kappa_0) - kappa_der * h_max
       soilpar(i)%capa_d =  kappa_0

       ! conductivity
       !h_max = h_max * 4.

       !h_max = -0.0001
       soilpar(i)%h_max2 = h_max   

       ! conductivity in the left and right
       kappa_0 = conduct(h_right, i, 0.)
       kappa_M = conduct(h_max, i, 0.)

       ! derivative of the condutivity in the left
       hh = abs( (h_right - h_left) / 200.) 
       kappa_der = (conduct(h_max+hh, i, 0.) - conduct(h_max-hh, i, 0.) ) / (2*hh)

       soilpar(i)%cond_a = -2.*(kappa_M - kappa_0) + kappa_der * h_max
       soilpar(i)%cond_b =  3.*(kappa_M - kappa_0) - kappa_der * h_max
       soilpar(i)%cond_d =  kappa_0

       ! correction moving the inflex point to the alpha-th part of (0, h_max)
       ! not working properly !!!
       alpha = 1.0

       A = (3*soilpar(i)%cond_a * alpha + soilpar(i)%cond_b ) !/h_max /h_max
       
       !B = ((1-alpha)**2 - 4 * alpha *(1-alpha) + alpha*alpha) ! / h_max /h_max
       B =  1 - 2 * alpha *(1-alpha)  ! / h_max /h_max
       
       !!soilpar(i)%cond_e =  - A / B

       soilpar(i)%cond_e =  0.
       
       !soilpar(i)%cond_e =  -1E-2 * kappa_0 / (h_max**4)
       
       ! write(*,'(a8, i2, 20es12.4)') 'cubic:',i, h_max, val_max, hh, kappa_der, kappa_0, kappa_M, &
       !      soilpar(i)%cond_a, soilpar(i)%cond_b, soilpar(i)%cond_d,soilpar(i)%cond_e

       ! write(30+i, *) h_right, kappa_0
       ! write(30+i, *) h_max, kappa_M

       ! write(40+i, *) i, h_max, val_max, kappa_M, '   jde9393h3h4343'

       ! write(50+i, *) h_max-hh, kappa_M - kappa_der*hh
       ! write(50+i, *) h_max, kappa_M
       ! write(50+i, *) h_max+hh, kappa_M + kappa_der*hh
       ! write(50+i,*) 
    end do ! i 

    cubic_interpol = .true.
    print*,' # Porous media flow problem, use of the cubic interpolation!!'
    
  end subroutine porous_coeffs_cubic_interpol
  
  
  !> \brief Tracy relation \f[ \theta = f(pressure) \f]
  !!  \f[ \theta = \frac{1}{(1+(\alpha*h)^n)^m} \f]
  !!
  !<
  function gardner_wc(u, layer, z) result(theta)
    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting water content
    real :: theta
   
    real :: h
    
    
    h = u - z


    if (h >=0.0) then
        theta = soilpar(layer)%Ths
        RETURN
    else
        theta = soilpar(layer)%Thr+(soilpar(layer)%Ths - soilpar(layer)%Thr)*exp(soilpar(layer)%alpha*h)
    end if

  end function gardner_wc
  
  
  !> \brief so-called retention water capacity, it is a derivative to retention curve function
  !! \f[ E(h) = \left\{ \begin{array}{l l} \alpha (\theta_s - \theta_r) exp(\alpha h) ,  & \quad \mbox{$\forall$ $h \in (-\infty, 0 )$}\\ 0, & \quad \mbox{$\forall$ $h \in \langle 0, + \infty )$}\\ \end{array} \right. \f]
  !<
  function gardner_cap(u, layer, z) result(E)
    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting system capacity (elasticity)
    real :: E

    real :: a,  tr, ts, h    
        
    h = u - z

    if (h < 0) then
      a = soilpar(layer)%alpha
      tr = soilpar(layer)%Thr
      ts = soilpar(layer)%Ths
      E = a*(-tr + ts)*exp(a*h)
    else
      E = 0
      RETURN
    end if


  end function gardner_cap

  

  
  !> \brief Tracy's fucntion for unsaturated hydraulic conductivity
  !! \f[   K(h) = K_s exp(\alpha h),  &  \mbox{$\forall$  $h \in$ $(-\infty,0)$}\\ K_s, & \mbox{$\forall$   $h \in$ $\langle 0, +\infty)$}\\ \end{array} \right. \f]
  !<
  function gardner_cond(u, layer, z) result(K)
    implicit none
    !> solution
    real, intent(in) :: u
    !> material id
    integer, intent(in) :: layer
    !> geodetic head
    real, intent(in) :: z
    !> resulting hydraulic conductivity
    real :: K
    real :: h
  
    h = u - z
    
    if (h >= 0) then
      K = soilpar(layer)%Ks
    else
      K = soilpar(layer)%Ks*exp(soilpar(layer)%alpha*h)
    end if

  end function gardner_cond
  
end module porous_fnc

!> Tracy (2006) analytical model (use Tracy's functions from porous_fnc
module re_analytical

  public :: tracy_fnc
  private :: tracy_bc

  contains

   !> Boundary condition for Tracy's analytical solution
    subroutine tracy_bc(bcval, x_coord, width, hinit)
      use porous_fnc
      implicit none

      !> output boundary value
      real, intent(out) :: bcval
      !> with - the width of the domain, coord - the x coordinate at the boundary, hihit - initial condition (should be very dry))
      real, intent(in) :: width, x_coord, hinit

      real :: hbar, alpha
    
      alpha = soilpar(1)%alpha
    
      hbar = 1 - exp(alpha*hinit)

      bcval = 1.0/alpha*log(exp(alpha*hinit) + hbar*sin(4*atan(1.0)*x_coord/width))


    end subroutine tracy_bc
  
   !> analytical solution to the transient 2D Richards' equation based on (Tracy, 2006)
    subroutine tracy(hinit, coord, t, width, length, h)
      use porous_fnc
      implicit none

      !> initial state (must be constant)
      real, intent(in) :: hinit
      !> point coordinates
      real, dimension(:), intent(in) :: coord
      !> simulation time
      real, intent(in)               :: t, width, length
      !> solution
      real, intent(out)              :: h



!--------local variables--------------------
      real :: lambda
      real :: c
      real :: gamma
      real :: phi
      real :: hbar
      real :: beta
      real :: ho
      real :: hr
      real :: suma
      real :: hss
      real :: tmp
      real :: alpha
      real :: a
      real :: L, absval
      integer :: i
      
      if (abs(t) < epsilon(t)) then
        if (abs(coord(2)-length) < epsilon(length)) then
          call tracy_bc(h, coord(1), width, hinit)
        else
          h=hinit
        end if
        RETURN
      end if
  

      a = width

      L = length


      alpha = soilpar(1)%alpha

      ho = 1-exp(alpha*hinit)

      beta = sqrt(alpha**2/4 + (4*atan(1.0)/a)**2)

      c = alpha*(soilpar(1)%ths-soilpar(1)%thr)/soilpar(1)%Ks

      hss = ho*sin(4*atan(1.0)*coord(1)/a)*exp(alpha/2*(L-coord(2)))*sinh(beta*coord(2))/sinh(beta*L)

      suma = 0

      i = 0

      !print*,'_______________________'
      do
        if (i<huge(i)/10) then
          i = i+1
        end if
        tmp = suma
        lambda = i*4*atan(1.0)/L
        gamma = 1/c*(beta*beta + lambda*lambda)
        tmp = ((-1)**i)*lambda/gamma*sin(lambda*coord(2))*exp(-gamma*t)
        if (i==1) absval=abs(1/c*lambda/gamma) * abs(exp(-gamma*t))
        suma = suma + tmp

        !write(*,'(a8, 60es12.4)') 'suma = ',t, suma,tmp, abs(1/c*lambda/gamma) * abs(exp(-gamma*t)), absval*epsilon(tmp)
        
        if (abs(1/c*lambda/gamma) * abs(exp(-gamma*t)) < absval*epsilon(tmp) ) then
          EXIT
        end if
      
      end do

      !print*,'finished = ',suma,tmp, abs(1/c*lambda/gamma) * abs(exp(-gamma*t)), absval*epsilon(tmp) 
        

      phi = 2*ho/(L*c)*sin(4*atan(1.0)*coord(1)/a)*exp(alpha/2*(L-coord(2)))*suma


      hbar = phi + hss

      ! print*
      ! print*,'xi, t = ', coord(1:2), t
      ! print*, 'hinit =', hinit
      ! print*,' width, length,  = ', width, length
      ! print*,'Ks = ', soilpar(1)%Ks
      ! print*,'Thr = ', soilpar(1)%thr
      ! print*,'Ths = ', soilpar(1)%ths
      ! print*, 'alpha = ', alpha, abs(1/c*lambda/gamma) * abs(exp(-gamma*t))
      ! print*, 'hinit = ', hinit
      ! print*, 'hbar = ', hbar, phi, hss, suma
      ! print*,'::A', (alpha*hinit)
      ! print*,'::B', (exp(alpha*hinit)+hbar)
      ! print*,'::C', log(exp(alpha*hinit)+hbar)

      
      h = 1/alpha*log(exp(alpha*hinit)+hbar)

      
      ! print*,'exact', coord(1:2), t, h
      ! write(23, *)  coord(1:2),  h, t
      ! print*,'______________________________'

  end subroutine tracy

end module re_analytical


!> definition of reconstruction  used in the model of porous media flow
module porous_data_module
  use paramets
  
  implicit none

  type :: porous_data_str
     integer :: npoin
     real, dimension(:, :), allocatable :: wR     ! node values of the reconstructed solution
     real, dimension(:, :), allocatable :: mater   ! node values of the material parameters
  end type porous_data_str

  ! characteristic parameters of the flow
  type :: porous_data_evals
     real :: water_content            ! at t_{m}^-
     real :: water_content_old        ! at t_{m-1}^-
     real :: water_content_back       ! at t_{m-1}^+
     real :: water_content_init       ! at t_{0}^-
     real :: water_content_losses     ! 
     real :: DBC_inaccuracy
     real :: DBC_inaccuracy_Time
     real :: flow_inlet
     real :: flow_bound
     real :: flow_inlet_Time
     real :: flow_bound_Time
     real :: flow_bilance
     real :: algeb_resid
     real :: algeb_resid_tot
  end type porous_data_evals

  type( porous_data_str) :: porous_data
  type( porous_data_evals) :: porous_evals

  contains

    ! initizalization of  characteristic parameters of the flow
    subroutine porous_data_evals_init( )

      ! DONE directly in  Eval_Characteristic_paramets
      !porous_evals%water_content = 0.
      !porous_evals%water_content_init = 0.
      !porous_evals%flow_inlet = 0.
      !porous_evals%flow_bound = 0.
      !porous_evals%flow_inlet_Time = 0.
      !porous_evals%flow_bound_Time = 0.
      !porous_evals%flow_inlet = 0.
      !porous_evals%flow_bilance = 0.
    end subroutine porous_data_evals_init
    
      
    !> initialization
    subroutine porous_data_str_init(npoin, ndimL)
      integer, intent(in) :: npoin, ndimL

      !print*,'______________________________________________________________'
      !print*
      !print*,'###########   porous_data_str_init', npoin, ndimL
      !print*,'______________________________________________________________'

      porous_data%npoin = npoin

      if(allocated(porous_data%wR)) then
         if(size (porous_data%wR, 1) /= npoin) then
            deallocate(porous_data%wR)
            allocate(porous_data%wR(1:npoin, 0:ndimL) , source = 0.0  )

            deallocate(porous_data%mater)
            allocate(porous_data%mater(1:npoin, 0:iRe), source = 0.0  )
         endif
      else
         allocate(porous_data%wR(1:npoin, 0:ndimL), source = 0.0  )
         allocate(porous_data%mater(1:npoin, 0:iRe), source = 0.0  )
      endif


    end subroutine porous_data_str_init

    !> initialization
    subroutine porous_data_str_deallocate( )
      !print*,'______________________________________________________________'
      !print*,'###########   porous_data_str_deallocate'
      !print*,'______________________________________________________________'

      deallocate(porous_data%wR )
      deallocate(porous_data%mater )

    end subroutine porous_data_str_deallocate

end module porous_data_module

!> definition of model of porous media flow
module porous_mod

  ! use main_data
  use model_mod
  !  use f_mapping
  !  use mesh_oper
  !  use define_state
  !  use blocks_integ
  !  use model3DNS
  !  use model2DNS
  !  use modelTurb2e
  !  use modelFE
  !  use modelLaplace
  use porous_fnc

  implicit none

  type, EXTENDS( Model_t ), public :: Porous_t
     !   real :: kappa, kappa1
     !   real :: Pr


contains
  procedure :: init => initPorous

end type Porous_t


contains
  !> initialization of the NS case
  !> isca, t2 not used
subroutine initPorous( this, Re, isca, t2)
 class (Porous_t), intent(inout) :: this
 real, intent(in) ::Re
 integer, intent(in), optional :: isca
 real, intent(in), optional :: t2

 !stop 'initPorous not implemented yet'
 !??? ndim?
 this%ndim = 1 ! is not sufficient for the global setting, set in readModelData !!
 this%convective = .false.

 this%subdomainRHS = .false.
 this%rhsTestFunDerivative = 0 ! rhs: \int f*phi dx
 ! is the exact solution a priori known?
 this%known_sol = .false.

 this%Re = Re
 this%icase = isca
 this%param1 = t2

 this%ireac = 0

 if( Re > 0.) then
    this%Re1  = 1./ this%Re
 else
    this%Re1 = 0.
 endif

 ! Tracy
 if (this%icase== 6)  USE_TRACY=.true.

 
 ! initialization of the data coefficients
 call init_porous_coeffs()

 ! cubic interpolation of coefficients near h<0
 call porous_coeffs_cubic_interpol(this%param1)
 print*,' # Porous media flow problem cubic remedy of capacity and conductivity'

 ! numerical integration of \f$ \theta(h) \f$ for the water content
 call porous_water_content_integ(this%icase )
 
 ! spline reconstruction for numerical integration of \f$ \theta(h) \f$ for the water content
 call porous_water_content_spline(  )
 print*,' # Porous media flow problem: cubic interpolation of the water content'

 !has to be modified !!!
 select case (this%icase)
 case(1)   ! damp (HRAZ),  linear test case
    this%idiff = 2
    this%iconv = 0
    this%iexact = 2
    this%conv_rat = 1.0

    !this%varying_time_term = .true. !  here only for the test

 case(2)   ! Forcheimer  damp (HRAZ),  Forchheimer 2-term law, only test case
    this%idiff = 2
    this%iconv = 0
    this%iexact = 2
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in fron of the time derivative

 case(3)   ! unsaturated flow through damp (HRAZ),  REAL data by M. Kuraz
    this%idiff = 3
    this%iconv = 0
    this%iexact = 2  
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in fron of the time derivative

 case(4)   ! test nonlinear case:  2u u_t - (u^2 u_x)_x = 0
    this%idiff = 4
    this%iconv = 0
    this%iexact = 3
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in fron of the time derivative

 case(5)   ! unsaturated flow 1D Wett
    this%idiff = 5
    this%iconv = 0
    this%iexact = 5  
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in fron of the time derivative


 case(6)   ! Tracy benchmark
    this%idiff = 6
    this%iconv = 0
    this%iexact = 6  
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in fron of the time derivative

 case(7)   ! test case for emission BC (seepgae BC)
    this%idiff = 5
    this%iconv = 0
    this%iexact = 7
    this%conv_rat = 1.0
    
    this%varying_time_term = .true. ! term in front of the time derivative

 case(8)   ! valcovazk
    this%idiff = 7
    this%iconv = 0
    this%iexact = 8  
    this%conv_rat = 1.0

    this%varying_time_term = .true. ! term in front of the time derivative


 case(9:71)
    print*,' UNKNOWN TYPE of scalar%icase !!!'
    print*,'Add definitions in:'
    print*,'                    o_porous.f90:    subroutine initPorous '
    print*,'                    modelPorous.f90: function Eval_Diff_Porous_Coeffs'
    print*,'                    modelPorous.f90: subroutine Set_Time_Matrix_porous'
    print*,'                    modelPorous.f90: subroutine Exact_Porous'
    print*,'                    problem.f90:      subroutine Setting_of_space_variable_coeffs'
    stop

 case(72)   ! porous media flow,  Forchheimer 2-term law
    this%idiff = 13
    this%iconv = 0
    this%iexact = 65
    this%conv_rat = 1.0
 case(73:)
    print*,' UNKNOWN TYPE of this%isca !!!'
    stop

 end select



 !    this%kappa = 1.4
 !    this%kappa1 = this%kappa - 1.
 !    this%Pr = 0.72   ! Prandtl number (almost all gases)
 !
 !
 !    if ( Re == 0.) then
 !        print*,' # Compressible Euler equations'
 !        this%Re1 = 0.
 !     elseif  ( Re > 0.) then
 !        this%Re1 = 1./this%Re
 !        print*,' # Compressible Navier-Stokes equations, Re=',this%Re
 !
 !     else
 !        print*,'# Reynolds number is negative',this%Re,' STABILIZATION !!!'
 !        !stop
 !     endif

end subroutine initPorous

end module porous_mod



