! 
! Copyright (c) 2008 Kristopher L. Kuhlman (kuhlman at hwr dot arizona dot edu)
! 
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
! 
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
! 
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! 
! Malama, B., A. Revil, and K.L. Kuhlman, 2009. A semi-analytical solution 
!  for transient streaming potentials associated with confined aquifer 
!  pumping tests, Geophysical Journal International, 176(3), 1007–1016.
! http://dx.doi.org/10.1111/j.1365-246X.2008.04014.x
!
! Malama, B., K.L. Kuhlman, and A. Revil, 2009. Theory of transient 
!  streaming potentials associated with axial-symmetric flow in unconfined 
!  aquifers, Geophysical Journal International, 179(2), 990–1003.
! http://dx.doi.org/10.1111/j.1365-246X.2009.04336.x
!

! $Id: laplace_hankel_solutions.f90,v 1.3 2009/03/15 03:38:14 kris Exp kris $

module three_layer
  implicit none

  private
  public :: unconfined_aquifer_sp

contains

  function unconfined_aquifer_sp(dum) result(fp)
    use shared_data, only :  stehfestN,tsval,zd,rd,bd1,bd3,&
         & layer,sigmaD,kappa,alphaD
    use constants, only : EP,DP,EONE, EZERO
    use invlap_fcns, only : stehfest_setup, stehfest

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none
    real(DP), intent(in) :: dum  ! scalar integration variable (Hankel parameter)
    real(DP) :: fp
    integer :: np
    real(EP), dimension(1:stehfestN) :: p,fpep
    real(EP), dimension(0:1,1:stehfestN) :: neuman
    real(EP), dimension(2) :: g,f,A
    real(EP) :: delta
    real(EP), dimension(1:2,1:stehfestN) :: xi
    real(EP), dimension(1:5) :: coshx,sinhx

    np = stehfestN
    call stehfest_setup(tsval,p(1:np))

    coshx(1) = cosh(real(dum*bd1,EP)); sinhx(1) = sinh(real(dum*bd1,EP))
    coshx(2) = cosh(real(dum*bd3,EP)); sinhx(2) = sinh(real(dum*bd3,EP))
    coshx(3) = cosh(real(dum,EP));     sinhx(3) = sinh(real(dum,EP))
    coshx(4) = cosh(real(dum*zd,EP));  sinhx(4) = sinh(real(dum*zd,EP))
    coshx(5) = cosh(real(dum*(1.0_DP - zd),EP))
    sinhx(5) = sinh(real(dum*(1.0_DP - zd),EP))    

    g(1) = coshx(1)*coshx(5) + real(sigmaD(1),EP)*sinhx(1)*sinhx(5)
    g(2) = coshx(2)*coshx(4) + real(sigmaD(3),EP)*sinhx(2)*sinhx(4)

    f(1) = (coshx(2)*coshx(3) + real(sigmaD(3),EP)*sinhx(2)*sinhx(3))/coshx(1)
    f(2) = sigmaD(1)*coshx(3)*sinhx(1) + sinhx(3)*coshx(1)

    A(1) = coshx(1)*coshx(2) + real(sigmaD(1)*sigmaD(3),EP)*sinhx(1)*sinhx(2)
    A(2) = real(sigmaD(1),EP)*sinhx(1)*coshx(2) + real(sigmaD(3),EP)*coshx(1)*sinhx(2)

    delta = abs(A(1)*sinhx(3) + A(2)*coshx(3))

    neuman(0,1:np) = neuman72(p(:),dum,0.0_DP,kappa,alphaD)
    neuman(1,1:np) = neuman72(p(:),dum,1.0_DP,kappa,alphaD)

    xi(1,1:np) = real(sigmaD(3),EP)*neuman(0,:)*sinhx(2)
    xi(2,1:np) = neuman(1,:)*(real(sigmaD(1),EP)*sinhx(1) - p(:)/(dum*alphaD)*coshx(1))

    if(layer(1)) then  ! top unpumped vadose zone
       fpep(1:np) = ((xi(1,:) + xi(2,:)*f(1))/delta - neuman(1,:)/coshx(1))* &
            & cosh(real(dum*(1.0_DP + bd1 - zd),EP))
    elseif(layer(2)) then  ! middle pumped unconfined layer
       fpep(1:np) = (xi(1,:)*g(1) + xi(2,:)*g(2))/delta - &
            & neuman72(p(:),dum,zd,kappa,alphaD)
    else ! lower no-flow aquiclude
       fpep(1:np) = (xi(2,:) - neuman(0,:)*f(2))/delta*cosh(real(dum*(bd3 + zd),EP))
    end if

    fp = dum*dbesj0(dum*rd)*stehfest(fpep(1:np),tsval)


  end function unconfined_aquifer_sp

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! this function is modified from the leaky-unconfined code
  !! it is only designed to be called from the one_aquifer_sp routine above
  pure function neuman72(p,dum,zd,kappa,alphaD) result(epfp)
    use constants, only : EP,DP,RTWO,ETWO,PIEP,EONE

    implicit none
    real(DP), intent(in) :: dum,zd,kappa,alphaD
    real(EP), intent(in), dimension(:) :: p
    real(EP), dimension(size(p)) :: theis,epfp,eta,xi
    integer :: np

    np = size(p)

    eta(1:np) = sqrt((p(:) + real(dum**2,EP))/real(kappa,EP))
    xi(1:np) = eta(:)*real(alphaD,EP)/p(:)
    theis(1:np) = ETWO/(p(:)*(p(:) + real(dum**2,EP)))

    epfp(1:np) = theis(:)*(EONE - cosh(eta(:)*real(zd,EP))/ &
         & (cosh(eta(:)) + xi(:)*sinh(eta(:))))

  end function neuman72

end module three_layer