! 
! Copyright (c) 2008 Kristopher L. Kuhlman (klkuhlm at sandia dot gov)
! 
! 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., K.L. Kuhlman, and W. Barrash, 2007. Semi-analytical solution
!  for flow in leaky unconfined aquifer-aquitard systems, Journal of
!  Hydrology, 346(1-2), 59–68.
! http://dx.doi.org/10.1016/j.jhydrol.2007.08.018
!
! Malama, B., K.L. Kuhlman, and W. Barrash, 2008. Semi-analytical solution
!  for flow in a leaky unconfined aquifer toward a partially penetrating
!  pumping well, Journal of Hydrology, 356(1-2), 234–244.
! http://dx.doi.org/10.1016/j.jhydrol.2008.03.029
!
	
! $Id: utility.f90,v 1.5 2007/02/19 02:21:26 kris Exp kris $

!! this file has _four_ modules in it 

!! ######################################################################
!! ######################################################################

module constants

  ! integers with range of 9 and 18 orders of magnitude
  integer, parameter :: I4B = selected_int_kind(r=9)   ! std 32-bit integer

  ! real with range 300 orders of mag, 15 sig figs (8 on both g95 & ifort)
  integer, parameter :: DP = selected_real_kind (p=15,r=300)

  !! extended range internal variables (10 on g95, 16 on ifort)
  integer, parameter :: EP = selected_real_kind(r=3000)

  ! useful? constants related to pi and ln
  ! calculated to precision=33 using ifort or mathematica
  real(DP), parameter :: PI =        3.141592653589793238462643383279503_DP
  real(DP), parameter :: TWOPI =     6.28318530717958647692528676655901_DP
  real(DP), parameter :: PIOV2 =     1.57079632679489661923132169163975_DP
  real(DP), parameter :: PIOV3 =     1.04719755119659774615421446109317_DP
  real(DP), parameter :: PIOV4 =     0.785398163397448309615660845819876_DP
  real(DP), parameter :: ONEOVPI =   0.318309886183790671537767526745029_DP
  real(DP), parameter :: LN10  =     2.30258509299404568401799145468436_DP
  real(DP), parameter :: LN2 =       0.693147180559945309417232121458177_DP
  real(DP), parameter :: LNPIOV2 =   0.451582705289454864726195229894882_DP 
  real(DP), parameter :: EULER =     0.5772156649015328606065120900824025_DP
  real(DP), parameter :: PISQ =      9.86960440108935861883449099987615_DP
  real(DP), parameter :: SQRTPI =    1.77245385090551602729816748334115_DP
  real(DP), parameter :: SQRTPIOV2 = 1.25331413731550025120788264240552_DP
  real(DP), parameter :: INVLOG10 =  0.434294481903251827651128918916605_DP

  real(EP), parameter :: PIEP = 3.141592653589793238462643383279503_EP

  ! a very small number (if < SMALL, it is effectively zero, compared to one)
  real(DP), parameter :: SMALL = epsilon(1.0_DP) ! approx 2.22E-16

  ! the largest number which can be represented accurately
  real(DP), parameter :: LARGE = huge(1.0_DP)  ! approx 1.79E+308

  ! commonly(?) used constants
  real(EP), parameter :: EONE = 1.0_EP
  real(EP), parameter :: EZERO = 0.0_EP
  real(EP), parameter :: EHALF = 0.5_EP
  real(EP), parameter :: ETWO = 2.0_EP
  complex(DP), parameter :: CZERO = (0.0_DP, 0.0_DP)
  complex(DP), parameter :: CONE  = (1.0_DP, 0.0_DP)
  complex(DP), parameter :: CTWO =  (2.0_DP, 0.0_DP)
  complex(DP), parameter :: EYE = (0.0_DP, 1.0_DP)
  real(DP),  parameter :: RZERO = 0.0_DP
  real(DP),  parameter :: RONE  = 1.0_DP
  real(DP),  parameter :: RTWO =  2.0_DP
  real(DP),  parameter :: HALF =  0.5_DP

end module constants

!! ######################################################################
!! ######################################################################

module shared_data
  use constants, only : DP,EP
  implicit none

  integer, parameter :: MAXNLAYERS = 3

  ! variables to pass using modules

  ! parameters to change between calculations
  real(DP), save :: tsval   ! time - passed through hankel integral to Lap xform

  ! aquifer parameters (primary - read from input)
  real(DP), save, dimension(MAXNLAYERS) :: b,Kr,Kz,Ss
  real(DP), save :: Sy,zd,rd, lambda, ld,dd

  ! aquifer parameters (secondary - computed from above)
  real(DP), save :: kappa,alphaDy,bd
  real(DP), save, dimension(MAXNLAYERS) :: alphar,alphaz,alphaDr,alphaDz,KDz

  real(DP), save :: ztop, zbot, bd2,bd3
  logical, save :: piezometer
  logical, dimension(MAXNLAYERS) :: layer

  ! inverse Laplace transform parameters
  real(DP) :: lapalpha, laptol
  integer :: lapM

  ! Gauss Quadrature local abscissa, weights, and abscissa in global coords
  real(EP), allocatable :: obsGQx(:), obsGQw(:), obsGQy(:)
  complex(EP), allocatable :: obsGQz(:,:)

end module shared_data

!! ######################################################################
!! ######################################################################

module complex_hyp_trig
  use constants, only : EP

  implicit none

  private
  public :: cecosh, cesinh, cecoth, cesinhcosh

contains

  !!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  !! extended-precision complex cosh
  pure elemental function cecosh(z) result(c)
    use constants, only : EHALF
    complex(EP), intent(in) :: z
    complex(EP) :: c

    c = EHALF*(exp(z) + exp(-z))

  end function cecosh

  !! extended-precision complex sinh
  pure elemental function cesinh(z) result(s)
    use constants, only : EHALF
    complex(EP), intent(in) :: z
    complex(EP) :: s

    s = EHALF*(exp(z) - exp(-z))

  end function cesinh

  !! extended-precision complex sinh and cosh 
  !! for same argument
  pure elemental subroutine cesinhcosh(z,s,c)
    use constants, only : EHALF
    complex(EP), intent(in) :: z
    complex(EP), intent(out) :: s,c
    complex(EP) :: expplus, expmin
    
    expplus = exp(z); expmin = exp(-z)
    c = EHALF*(expplus + expmin)
    s = EHALF*(expplus - expmin)

  end subroutine cesinhcosh

  !! extended-precision complex coth
  pure elemental function cecoth(z) result(c)
    use constants, only : EONE,ETWO
    complex(EP), intent(in) :: z
    complex(EP) :: c, temp

    temp = exp(-ETWO*z)
    c = (EONE + temp)/(EONE - temp)

  end function cecoth

end module complex_hyp_trig

module invlap_fcns
implicit none
private
public :: invlap_setup, check_nan

contains

  !! consolidate calculation of laplace parameter vector
  subroutine invlap_setup(p,tee)
    use constants, only : DP,EP,PIEP,ETWO
    use shared_data, only : lapM,lapAlpha,lapTol,tsval
    
    complex(EP), intent(out), dimension(2*lapM+1) :: p
    real(DP), intent(out) :: tee
    integer :: np,i
    real(EP), dimension(0:2*lapM) :: run

    run = real((/ (i,i=0,2*lapM) /),EP)
    np = 2*lapM+1
    tee = 2.0_DP*tsval

    p(1:np) = cmplx(real(lapalpha,EP) - log(real(laptol,EP))/(ETWO*tee), run*PIEP/tee, EP)

    ! check for potential problems
    if(minval(abs(p)) < lapalpha) then
       stop 'decrease alpha parameter for de Hoog routine'
    end if
    
  end subroutine invlap_setup

  subroutine check_nan(fp,np,M)
    use constants, only : EP
    
    complex(EP), intent(in), dimension(:) :: fp
    integer, intent(out) :: np, M
    logical, dimension(size(fp)) :: mask
    integer :: nfp, j

    nfp = size(fp)
    mask = .false.
    np = nfp
    M = (np-1)/2

    where(fp /= fp)
       mask = .true.
    end where
    
    loop: do j=1,nfp
       if(mask(j)) then
          np = j
          exit loop
       end if
    end do loop

    if(np /= nfp) then
       if(np <= 4) then
          print *, 'insufficient terms for useful Laplace inversion',np
       else
          print *, 'reduced np:',nfp,'->',np
       end if
       
       if(mod(np,2) == 1) then
          !! if new np is odd
          M = (np-1)/2
       else
          !! if new np is even, lose additional term
          M = (np-2)/2
          np = np-1
       end if
    end if

    
  end subroutine check_nan
  
end module invlap_fcns



