! 
! 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., 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: laplace_hankel_solutions.f90,v 1.16 2008/03/10 21:47:14 kris Exp kris $

!! ######################################################################
!! there is a module in this file for each Hankel-Laplace space solution
!! each solution is vectorized with respect to the Laplace parameter
!! ######################################################################

module confined_partialpen
  implicit none

  private
  public :: partialpen_laphank, partialpen_noleak

  !! change here if you want to use a different one of the
  !! three (hopefully equivalent) versions of the Hantush solution

  interface partialpen_laphank
     module procedure partialpen_laphank1
  end interface
  
contains

  !! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !! this is the formulation of ud given in the text 
  !! it is the most mathematically simplified form
  function partialpen_laphank1(eta1,zd,p) result(fp)
    use constants, only : DP,EP,ETWO,EONE,EZERO
    use shared_data, only :ld,dd,alphaDz
    use complex_hyp_trig, only : cecosh,cesinh

    implicit none

    complex(EP), intent(in), dimension(:) :: eta1,p
    real(EP), intent(in) :: zd
    real(EP) :: epdd,epld, zeta
    complex(EP), dimension(size(p)) ::fp, deltaud
    integer :: np
#ifdef DEBUG
    integer :: k
#endif

    np = size(p); epdd = real(dd,EP); epld = real(ld,EP)

    deltaud(1:np) = (cesinh(eta1(:)*epdd)*cecosh(eta1(:)*(EONE + zd)) &
         & + cesinh(eta1(:)*(EONE - epld))*cecosh(eta1(:)*zd))/cesinh(eta1(:))

    if(EZERO >= zd .and. zd > -epdd) then  ! above pumping well screen 
       zeta = epdd + zd
    elseif(-epdd >= zd .and. zd >= -epld) then  ! in pumping well screen
       zeta = EZERO
    elseif(-epld > zd .and. zd >= -EONE) then  ! below pumping well screen
       zeta = epld + zd
    else
       stop 'PARTIALPEN_LAPHANK: invalid value of zd'
    end if

    fp(1:np) = ETWO*(cecosh(eta1(:)*zeta) - deltaud(1:np)) &
         & /(p(:)*eta1(:)*eta1(:)*alphaDz(1)*(epld-epdd))

#ifdef DEBUG
    open(unit=77,file='debug_partialpen.dat',status='replace',action='write')
    write(77,'(A,2(ES23.14))') '# deltaud, left',real(p(1)),aimag(p(2))

    do k=1,np
       write(77,777) k,real(deltaud(k)),aimag(deltaud(k)),real(cecosh(eta1(k)*zeta)), &
            & aimag(cecosh(eta1(k)*zeta)),real(p(k)*eta1(k)**2 *alphaDz(1)*(epld-epdd)),&
            & aimag(p(k)*eta1(k)**2 *alphaDz(1)*(epld-epdd)),real(fp(k)),aimag(fp(k))
    end do

777 format(I3,8ES17.8E2)
    close(77)
#endif

  end function partialpen_laphank1

  !! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !! this is the unsimplified form of ud given in the text 
  !! it is in terms of the 6 unknown coefficients which are solved for
  function partialpen_laphank2(eta1,zd,p) result(fp)
    use constants, only : DP,EP,ETWO,EONE,EZERO
    use shared_data, only :ld,dd,alphaDz
    use complex_hyp_trig, only : cesinh

    implicit none

    complex(EP), intent(in), dimension(:) :: eta1,p
    real(EP), intent(in) :: zd
    real(EP) :: epdd,epld
    complex(EP), dimension(size(p)) :: fp, denom, sinh1
    complex(EP), dimension(3,size(p)) :: A,B
    complex(EP), dimension(2,size(p)) :: wd
    integer :: np
#ifdef DEBUG
    integer :: k, zdi
    character(3) :: chnp,chz
    character(18) :: fmt
    complex(EP), dimension(size(p)) :: tmp
    zdi = nint(zd)
#endif

    np = size(p); epdd = real(dd,EP); epld = real(ld,EP)
    sinh1(1:np) = cesinh(eta1)

    wd(1,:) = (exp(eta1)*cesinh(eta1*epdd) + cesinh(eta1*(EONE - epld)))/sinh1(:)
    wd(2,:) = (exp(-eta1)*cesinh(eta1*epdd) + cesinh(eta1*(EONE - epld)))/sinh1(:)

#ifdef CONSISTENCY
    write(chz,'(I2.2)') abs(zdi)
    open(unit=77,file='consistency_pp'//trim(chz)//'.out',action='write',status='replace')
    write(chnp,'(I3.3)') 2*np
    fmt = '('//chnp//'ES24.14E4)'
    write(77,'(A)') 'all values should be zero'
#endif

    denom = p*eta1**2*alphaDz(1)*(epld - epdd)

    ! above pumping well screen 
    ! use A1 and B1
    A(1,:) = (exp(eta1*epdd) - wd(1,:))/denom
    B(1,:) = (exp(-eta1*epdd) - wd(2,:))/denom

    ! in pumping well screen
    ! use A3 and B3
    A(3,:) = -wd(1,:)/denom
    B(3,:) = -wd(2,:)/denom

    ! below pumping well screen
    ! use A2 and B2
    A(2,:) = (exp(eta1*epld) - wd(1,:))/denom
    B(2,:) = (exp(-eta1*epld) - wd(2,:))/denom

#ifdef CONSISTENCY
    tmp = A(1,:) - B(1,:)
    write(77,'(A,ES10.4)') 'BC at zD=0 (A-15) <=',maxval(abs(tmp))
    write(77,fmt) tmp 

    tmp = A(2,:)*exp(-eta1) - B(2,:)*exp(eta1)
    write(77,'(A,ES10.4)') 'BC at zD=-1 (A-16) <=',maxval(abs(tmp))
    write(77,fmt) tmp 

    tmp = (A(1,:) - A(3,:))*exp(-eta1*epdd) + (B(1,:) - B(3,:))*exp(eta1*dD) - ETWO/denom
    write(77,'(A,ES10.4)') 'head continuity at zD=-dD (A-17) <=',maxval(abs(tmp))
    write(77,fmt) tmp 

    tmp = (A(1,:) - A(3,:))*exp(-eta1*epdd) - (B(1,:) - B(3,:))*exp(eta1*dD)
    write(77,'(A,ES10.4)') 'flux continuity at zD=-dD (A-18) <=',maxval(abs(tmp))
    write(77,fmt) tmp 

    tmp = (A(2,:) - A(3,:))*exp(-eta1*epld) + (B(2,:) - B(3,:))*exp(eta1*lD) - ETWO/denom
    write(77,'(A,ES10.4)') 'head continuity at zD=-lD (A-19) <=',maxval(abs(tmp))
    write(77,fmt) tmp 

    tmp = (A(2,:) - A(3,:))*exp(-eta1*epld) - (B(2,:) - B(3,:))*exp(eta1*lD)
    write(77,'(A,ES10.4)') 'flux continuity at zD=-lD (A-19) <=',maxval(abs(tmp))
    write(77,fmt) tmp 
    close(77)
#endif

    if(EZERO >= zd .and. zd > -epdd) then 
       fp(1:np) = A(1,:)*exp(eta1*zd) + B(1,:)*exp(-eta1*zd)
    elseif(-epdd >= zd .and. zd >= -epld) then 
       fp(1:np) = A(3,:)*exp(eta1*zd) + B(3,:)*exp(-eta1*zd) + ETWO/denom
    elseif(-epld > zd .and. zd >= -EONE) then 
       fp(1:np) = A(2,:)*exp(eta1*zd) + B(2,:)*exp(-eta1*zd)
    end if

  end function partialpen_laphank2

  !! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  !! this formulation is in terms of "simplified" expressions 
  !! for ud in each interval, found using Mathematica via
  !! FullSimplify[TrigToExp[ud],Assumptions]

  function partialpen_laphank3(eta1,zd,p) result(fp)
    use constants, only : DP,EP,ETWO,EONE,EZERO,EHALF
    use shared_data, only :ld,dd,alphaDz
    use complex_hyp_trig, only : cecoth

    implicit none

    complex(EP), intent(in), dimension(:) :: eta1,p
    real(EP), intent(in) :: zd
    real(EP) :: epdd,epld
    complex(EP), dimension(size(p)) ::fp, denom
    integer :: np

    np = size(p); epdd = real(dd,EP); epld = real(ld,EP)
    denom(1:np) = p(:)*eta1(:)*eta1(:)*alphaDz(1)*(epld-epdd)

    if(EZERO >= zd .and. zd > -epdd) then
       ! above pumping well screen 
       fp(1:np) = (exp(-(epdd+epld+zd)*eta1)*(exp(epdd*eta1) - &
            & exp(epld*eta1))*(exp(ETWO*eta1) + exp((epdd + epld)*eta1))*(&
            & EONE + exp(ETWO*zd*eta1))*(cecoth(eta1) - EONE))/(-ETWO*denom)
    elseif(-epdd >= zd .and. zd >= -epld) then
       ! in pumping well screen
       fp(1:np) = (ETWO*(EONE - exp(eta1)/4.0_EP*(-exp(-(epld+zd-EONE)*eta1)*(&
            &  exp(ETWO*(epld - EONE)*eta1) - EONE)*(EONE + exp(ETWO*zd*eta1)) + &
            & exp(-(EONE+epdd+zd)*eta1)*(exp(ETWO*epdd*eta1) - EONE)*(EONE + &
            & exp(ETWO*(EONE + zd)*eta1)))*(cecoth(eta1) - EONE)))/denom
    elseif(-epld > zd .and. zd >= -EONE) then
       ! below pumping well screen
       fp(1:np) = (exp(-(epdd+epld+zd)*eta1)*(exp(epdd*eta1) - exp(epld*eta1) + &
            & exp((ETWO*epdd + epld)*eta1) - exp((epdd + ETWO*epld)*eta1))*(&
            & EONE + exp(ETWO*(EONE + zd)*eta1))*(cecoth(eta1) - EONE))/(-ETWO*denom)
    end if

  end function partialpen_laphank3

  !##################################################
  function partialpen_noleak(dum) result(f)
    use constants, only : EP,DP
    use shared_data, only : lapalpha,laptol,lapM,tsval,rd,zd,piezometer,&
         & obsGQy,obsGQw,obsGQz,alphaDz
    use invlap, only : dehoog_invlap_ep
    use invlap_fcns

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none
    real(DP), intent(in) :: dum
    real(DP) :: f, tee
    integer :: np, k, ord, np2, M2
    complex(EP), dimension(1:2*lapM+1) :: p
    complex(EP), dimension(1:2*lapM+1) :: epfp,eta1

    np = 2*lapM+1
    call invlap_setup(p(1:np),tee)

    ! alphaDr(1) == 1
    eta1(1:np) = sqrt((p(:) + real(dum,EP)**2)/alphaDz(1))

    if(piezometer) then 
       ! one value of zd
       epfp(1:np) = partialpen_laphank(eta1(:),real(zd,EP),p(:))
    else
       ord = size(obsGQy)
       ! integrate over zd for each value of p using Gauss Quadrature
       do k=1, ord
          obsGQz(k,1:np) = partialpen_laphank(eta1(:),obsGQy(k),p(:))
       end do
       epfp(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
    end if

    call check_nan(epfp,np2,M2)

    ! argument to Hankel integral
    f = dum*dbesj0(dum*rd)*real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,epfp(1:np2),M2),DP)

  end function partialpen_noleak
end module confined_partialpen

!! ######################################################################
!! ######################################################################
module two_aquifer
  implicit none

  private
  public :: twoAquifer_leak

contains

  function twoAquifer_leak(dum) result(f)
    use shared_data, only :  lapalpha,laptol,tsval,lapM,zd,rd,piezometer,bd2,bd3,&
         & layer,obsGQy,obsGQw,obsGQz,alphaDz,alphaDr,alphaDy,Kdz,ld,dd
    use constants, only : EP,DP,PI,RONE,RTWO,RZERO,PIEP,ETWO,EONE,EHALF,EZERO
    use invlap, only : dehoog_invlap_ep    ! invlap.f90
    use complex_hyp_trig, only : cesinhcosh, cesinh
    use invlap_fcns
    use confined_partialpen, only : partialpen_laphank

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none
    real(DP), intent(in) :: dum  ! scalar integration variable (Hankel parameter)
    real(DP) :: f, tee
    integer :: np, ord, k, np2, M2
    complex(EP), dimension(1:2*lapM+1) :: p,xi,delta,fpep,feta,denom
    complex(EP), dimension(2,1:2*lapM+1) :: theta,gamma
    complex(EP), dimension(-1:1,1:2*lapM+1) :: ud
    complex(EP), dimension(3,1:2*lapM+1) :: eta
    complex(EP), dimension(8,1:2*lapM+1) :: coshx,sinhx
    complex(EP), dimension(6,1:2*lapM+1) :: c
#ifdef DEBUG
    integer :: j
    character(3) :: chnp
    character(18) :: fmt    
    complex(EP), dimension(1:2*lapM+1) :: tmp
    complex(DP), dimension(6,1:2*lapM+1) :: cm, diff
#endif

    np = 2*lapM+1
    ord = size(obsGQy)
    call invlap_setup(p(1:np),tee)

    forall (k=1:3)
       eta(k,1:np) = sqrt((p(:) + real(alphaDr(k)*dum**2,EP))/alphaDz(k))
    end forall
    xi(1:np) = eta(1,:)*real(alphaDy,EP)/p(:)
    forall (k=1:2)
       gamma(k,1:np) = eta(k+1,:)*real(Kdz(k+1),EP)/eta(k,:)
    end forall

    theta(1,1:np) = eta(3,:)*(real(bd3,EP) - real(bd2,EP)) + eta(2,:)*(real(bd2,EP) - EONE)
    theta(2,1:np) = eta(3,:)*(real(bd3,EP) - real(bd2,EP)) - eta(2,:)*(real(bd2,EP) - EONE)

    call cesinhcosh(eta(1,:) + theta(1,:),sinhx(1,:),coshx(1,:))
    call cesinhcosh(eta(1,:) - theta(1,:),sinhx(2,:),coshx(2,:))
    call cesinhcosh(eta(1,:) + theta(2,:),sinhx(3,:),coshx(3,:))
    call cesinhcosh(eta(1,:) - theta(2,:),sinhx(4,:),coshx(4,:))
    call cesinhcosh(eta(1,:) - theta(2,:),sinhx(7,:),coshx(7,:))
    call cesinhcosh(eta(1,:) + theta(2,:),sinhx(8,:),coshx(8,:))

    delta(1:np) =((EONE + gamma(1,:))*(EONE + gamma(2,:))*(xi(:)*sinhx(1,:) + coshx(1,:)) &
         & + (EONE - gamma(1,:))*(EONE + gamma(2,:))*(xi(:)*sinhx(2,:) + coshx(2,:)) &
         & + (EONE - gamma(1,:))*(EONE - gamma(2,:))*(xi(:)*sinhx(3,:) + coshx(8,:)) & 
         & + (EONE + gamma(1,:))*(EONE - gamma(2,:))*(xi(:)*sinhx(4,:) + coshx(7,:)))/ETWO

    if(any(delta /= delta)) print *, 'NaN in Delta'

    call cesinhcosh(eta(1,:),sinhx(5,:),coshx(5,:))
    call cesinhcosh(eta(3,:)*(real(bd3,EP) - real(bd2,EP)),sinhx(6,:),coshx(6,:))

    !! simplify ud for specific values of zd
    denom(1:np) = sinhx(5,:)*p*eta(1,:)**2*alphaDz(1)*(real(ld,EP)-real(dd,EP))
    ud(0,1:np) = ETWO*(cesinh(eta(1,:)*(EONE - real(dd,EP))) - &
         & cesinh(eta(1,:)*(EONE - real(ld,DP))))/denom
    ud(-1,1:np) = ETWO*(cesinh(eta(1,:)*real(dd,EP)) - &
         &  cesinh(eta(1,:)*real(ld,DP)))/(-denom)

    feta(1:np) = ud(0,:) - ud(-1,:)*(xi(:)*sinhx(5,:) + coshx(5,:))

    c(3,1:np) = -exp(eta(2,:)*real(bd2,EP))*feta(:)*(coshx(6,:) + gamma(2,:)*sinhx(6,:))/delta(:)
    c(4,1:np) = -exp(-eta(2,:)*real(bd2,EP))*feta(:)*(coshx(6,:) - gamma(2,:)*sinhx(6,:))/delta(:)

#ifdef CONSISTENCY
    open(unit=88,file='consistency.out',action='write',status='replace')
    write(chnp,'(I3.3)') np
    fmt = '('//chnp//'ES24.14E4)'
    write(88,'(A)') 'all values should be zero'
#endif

#ifdef MATRIXCHECK
    cm(1:6,1:np) = numerical_matrix_twoaquifer(p(:),eta(1:3,:),xi(:),gamma(1:2,:),ud(-1:0,:))
    open(unit=89,file='matrix_check.out',action='write',status='replace')

    diff(3:4,1:np) = cmplx(real(c(3:4,:)),aimag(c(3:4,:)),DP) - cm(3:4,:)
    write(89,*) 'c3 max relative difference',maxval(abs(diff(3,:)/c(3,:)))
    write(89,*) diff(3,:)/c(3,:)
    write(89,*) 'c4 max relative difference',maxval(abs(diff(4,:)/c(4,:)))
    write(89,*) diff(4,:)/c(4,:)
#endif

    !! things specific to only one layer
    if(layer(1)) then  ! upper unconfined pumped aquifer
       c(1,1:np) = -exp(eta(1,:))/ETWO*(ud(-1,:) - c(3,:)*(EONE + gamma(1,:))*exp(-eta(2,:))&
            & - c(4,:)*(EONE - gamma(1,:))*exp(eta(2,:)))
       c(2,1:np) = -exp(-eta(1,:))/ETWO*(ud(-1,:) - c(3,:)*(EONE - gamma(1,:))*exp(-eta(2,:))&
            & - c(4,:)*(EONE + gamma(1,:))*exp(eta(2,:)))

#ifdef MATRIXCHECK
       diff(1:2,1:np) = cmplx(real(c(1:2,:)),aimag(c(1:2,:)),DP) - cm(1:2,:)
       write(89,*) 'c1 max relative difference',maxval(abs(diff(1,:)/c(1,:)))
       write(89,*) diff(1,:)/c(1,:)
       write(89,*) 'c2 max realtive difference',maxval(abs(diff(2,:)/c(2,:)))
       write(89,*) diff(2,:)/c(2,:)
#endif

#ifdef CONSISTENCY
       tmp = -c(1,:)*(xi(:) + EONE) + c(2,:)*(xi(:) - EONE) 
       write(88,'(A,ES10.4)') 'water table BC at zD=0 -simplified ud0- (A-30) <=',&
            & maxval(abs(tmp - ud(0,:)))
       write(88,fmt) tmp - ud(0,:)
       ud(0,:) = partialpen_laphank(eta(1,:),EZERO,p(:))
       write(88,'(A,ES10.4)') 'water table BC at zD=0 (A-30) -new-new ud0- <=',&
            & maxval(abs(tmp - ud(0,:)))
       write(88,fmt) tmp - ud(0,:)
       tmp = -c(1,:)*exp(-eta(1,:)) - c(2,:)*exp(eta(1,:)) &
            & + c(3,:)*exp(-eta(2,:)) + c(4,:)*exp(eta(2,:))
       write(88,'(A,ES10.4)') 'head continuity at zD=-1 (A-33) :simplified ud-1: <=',&
            & maxval(abs(tmp  - ud(-1,:)))
       write(88,fmt) tmp - ud(-1,:)
       ud(-1,:) = partialpen_laphank(eta(1,:),-EONE,p(:))
       write(88,'(A,ES10.4)') 'head continuity at zD=-1 (A-33) :new-new ud-1: <=',&
            & maxval(abs(tmp  - ud(-1,:)))
       write(88,fmt) tmp - ud(-1,:)
       tmp = c(1,:)*exp(-eta(1,:)) - c(2,:)*exp(eta(1,:)) &
            & - gamma(1,:)*c(3,:)*exp(-eta(2,:)) + gamma(1,:)*c(4,:)*exp(eta(2,:))  
       write(88,'(A,ES10.4)') 'flux continuity at zD=-1 (A-34) <=', maxval(abs(tmp))
       write(88,fmt) tmp
#endif

       if(piezometer) then
          ud(+1,:) = partialpen_laphank(eta(1,:),real(zd,EP),p(:))
          fpep(1:np) =  ud(1,:) + c(1,:)*exp(eta(1,:)*real(zd,EP)) + &
               & c(2,:)*exp(-eta(1,:)*real(zd,EP))
       else
          do k=1,ord
             ud(+1,:) = partialpen_laphank(eta(1,:),obsGQy(k),p(:))
             obsGQz(k,1:np) = ud(1,:) + &
                  & c(1,:)*exp(eta(1,:)*obsGQy(k)) + c(2,:)*exp(-eta(1,:)*obsGQy(k))
          end do
          fpep(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
       end if

    elseif(layer(2)) then  ! middle confined unpumped aquitard
       if(piezometer) then
          fpep(1:np) = c(3,:)*exp(eta(2,:)*real(zd,EP)) + c(4,:)*exp(-eta(2,:)*real(zd,EP))
       else
          do k=1,ord
             obsGQz(k,1:np) = c(3,:)*exp(eta(2,:)*obsGQy(k)) + c(4,:)*exp(-eta(2,:)*obsGQy(k))
          end do
          fpep(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
       end if

    else ! lower confined unpumped aquifer

       c(5,1:np) = exp(eta(3,:)*real(bd2,EP))/ETWO*&
            & (c(3,:)*(EONE + EONE/gamma(2,:))*exp(-eta(2,:)*real(bd2,EP)) &
            & + c(4,:)*(EONE - EONE/gamma(2,:))*exp(eta(2,:)*real(bd2,EP)))
       c(6,1:np) = exp(-eta(3,:)*real(bd2,EP))/ETWO*&
            & (c(3,:)*(EONE - EONE/gamma(2,:))*exp(-eta(2,:)*real(bd2,EP)) &
            & + c(4,:)*(EONE + EONE/gamma(2,:))*exp(eta(2,:)*real(bd2,EP)))

#ifdef MATRIXCHECK
       diff(5:6,1:np) = cmplx(real(c(5:6,:)),aimag(c(5:6,:)),DP) - cm(5:6,:)
       write(89,*) 'c5 max relative difference',maxval(abs(diff(5,:)/c(5,:)))
       write(89,*) diff(5,:)/c(5,:)
       write(89,*) 'c6 max relative difference',maxval(abs(diff(6,:)/c(6,:)))
       write(89,*) diff(6,:)/c(6,:)
#endif

#ifdef CONSISTENCY
       tmp = c(5,:)*exp(-eta(3,:)*real(bd3,EP)) - c(6,:)*exp(eta(3,:)*real(bd3,EP))
       write(88,'(A,ES10.4)') 'no-flow BC at zD=-bD,3 (A-36) <=',maxval(abs(tmp))
       write(88,fmt) tmp
       tmp = c(3,:)*exp(-eta(2,:)*real(bD2,EP)) + c(4,:)*exp(eta(2,:)*real(bd2,EP)) &
            & - c(5,:)*exp(-eta(3,:)*real(bd2,EP)) - c(6,:)*exp(eta(3,:)*real(bd2,EP))
       write(88,'(A,ES10.4)') 'head continuity at zD=-bD,2 (A-37) <=',maxval(abs(tmp))
       write(88,fmt) tmp
       tmp = c(3,:)*exp(-eta(2,:)*real(bD2,EP)) - c(4,:)*exp(eta(2,:)*real(bd2,EP)) &
            & - c(5,:)*gamma(2,:)*exp(-eta(3,:)*real(bd2,EP)) &
            & + c(6,:)*gamma(2,:)*exp(eta(3,:)*real(bd2,EP))
       write(88,'(A,ES10.4)') 'flux continuity at zD=-bD,2 (A-38) <=',maxval(abs(tmp))
       write(88,fmt) tmp
#endif  

       if(piezometer) then
          fpep(1:np) = c(5,:)*exp(eta(3,:)*real(zd,EP)) + c(6,:)*exp(-eta(3,:)*real(zd,EP))
       else
          do k=1,ord
             obsGQz(k,1:np) = c(5,:)*exp(eta(3,:)*obsGQy(k)) + c(6,:)*exp(-eta(3,:)*obsGQy(k))
          end do
          fpep(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
       end if
    end if

    call check_nan(fpep,np2,M2)
    f = dum*dbesj0(dum*rd)* &
         & real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,fpep(1:np2),M2),DP)

#ifdef DEBUG
    write(*,*) '%%%%%'
    write(*,*) 'Re(p)',real(p(1))
    write(*,*) 'Im(p)',aimag(p(2))
    write(*,*) 'a', dum
    if(piezometer) then
       write(*,*) 'zd', zd
    else
       write(*,*) 'obsGQy',obsGQy
    end if
    write(*,*) 'rd', rd 
    write(*,*) 'alphaDr(2:3)', alphaDr(2:3)
    write(*,*) 'alphaDz(1:3)', alphaDz(1:3)

    open(unit=55,file='debug_3layer_1.dat',status='replace',action='write')
    open(unit=66,file='debug_3layer_2.dat',status='replace',action='write')
    write(55,'(A,2(ES23.14))') '# p, eta{1,2,3}, xi, gamma{1,2}',real(p(1)),aimag(p(2))
    write(66,'(A,2(ES23.14))') '# p, theta{1,2}, delta, ud{-1,0}',real(p(1)),aimag(p(2))

    do k=1,np
       write(55,555) k,(real(eta(j,k)),aimag(eta(j,k)),j=1,3), &
            & real(xi(k)),aimag(xi(k)),(real(gamma(j,k)),aimag(gamma(j,k)),j=1,2)
       write(66,556) k,(real(theta(j,k)),aimag(theta(j,k)),j=1,2),&
            & real(delta(k)),aimag(delta(k)),(real(ud(j,k)),aimag(ud(j,k)),j=-1,0), &
            & real(fpep(k)),aimag(fpep(k))

    end do

555 format(I3,12(ES15.6E3))
556 format(I3,12(ES16.5E4))
    close(55); close(66)
    stop 'debugging stop'
#endif

#ifdef CONSISTENCY
    close(88)
    stop 'consistency check stop'
#endif

#ifdef MATRIXCHECK
    close(89)
    stop 'matrix check stop'
#endif

  end function twoAquifer_leak


#ifdef MATRIXCHECK
  !! ######################################################################
  !! this function is only DOUBLE PRECISION (since that is what LAPACK is)
  !! ######################################################################
  !! it avoids evaluating sinh,cosh,tanh, or doing any algebra  

  function numerical_matrix_twoaquifer(p,eta,xi,gamma,ud) result(nc)
    use shared_data, only :  bd2,bd3
    use constants, only : EP,DP
    use invlap_fcns

    implicit none
    interface  ! F95 interface to F77 LAPACK routine
       subroutine ZGESVX(FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, &
            & IPIVOT, EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, &
            & BERR, WORK, RWORK, INFO)

         character(LEN=1), intent(in) :: FACT, TRANS
         character(LEN=1), intent(inout) :: EQUED
         integer, intent(in) :: N, NRHS, LDA, LDAF, LDB, LDX
         complex(8), intent(inout), dimension(LDA,N) :: A
         complex(8), intent(inout), dimension(LDAF,N) :: AF
         complex(8), intent(inout), dimension(LDB,NRHS) :: B
         complex(8), intent(out), dimension(LDX,NRHS) :: X
         integer, intent(out) :: INFO
         integer, intent(inout), dimension(N) :: IPIVOT
         real(8), intent(out) :: RCOND
         real(8), intent(out), dimension(NRHS) :: FERR, BERR
         real(8), intent(inout), dimension(N) :: R, C
         real(8) , intent(inout), dimension(2*N) :: RWORK
         complex(8), intent(inout), dimension(2*N) :: WORK
       end subroutine ZGESVX
    end interface

    character(1) :: eqed
    integer :: np, j, k, ierr
    complex(EP), intent(in), dimension(:) :: p,xi
    complex(EP), intent(in), dimension(:,:) :: gamma, eta
    complex(EP), intent(in), dimension(-1:,:) :: ud
    complex(DP), dimension(6,size(p)) :: nc
    real(DP), dimension(size(p)) :: rcond,ferr,berr

    complex(DP), dimension(6,size(p)) :: b
    complex(DP), dimension(6,6,size(p)) :: A,AF

    ! things needed for LAPACK routine
    complex(DP), dimension(12) :: work
    real(DP), dimension(12) :: rwork
    integer, dimension(6) :: ipiv
    real(DP), dimension(6) :: lpR,lpC

    np = size(p)

    ! initialize
    A(1:6,1:6,1:np) = (0.0_DP, 0.0_DP)
    b(3:6,1:np) =     (0.0_DP, 0.0_DP)

    ! EQ A-30
    A(1,1,1:np) = -(xi(:) + 1.0_DP)
    A(1,2,1:np) =  (xi(:) - 1.0_DP)
    b(1,1:np) = ud(0,:)

    ! EQ A-33
    A(2,1,1:np) = -exp(-eta(1,:))
    A(2,2,1:np) = -exp( eta(1,:))
    A(2,3,1:np) =  exp(-eta(2,:))
    A(2,4,1:np) =  exp( eta(2,:))
    b(2,1:np) = ud(-1,:)

    ! EQ A-34
    A(3,1,1:np) =  exp(-eta(1,:))
    A(3,2,1:np) = -exp( eta(1,:))
    A(3,3,1:np) = -exp(-eta(2,:))*gamma(1,:)
    A(3,4,1:np) =  exp( eta(2,:))*gamma(1,:)   

    ! EQ A-36
    A(4,5,1:np) =  exp(-eta(3,:)*bd3)
    A(4,6,1:np) = -exp( eta(3,:)*bd3)

    ! EQ A-37
    A(5,3,1:np) =  exp(-eta(2,:)*bd2)
    A(5,4,1:np) =  exp( eta(2,:)*bd2)
    A(5,5,1:np) = -exp(-eta(3,:)*bd2)
    A(5,6,1:np) = -exp( eta(3,:)*bd2)

    ! EQ A-38
    A(6,3,1:np) =  exp(-eta(2,:)*bd2)
    A(6,4,1:np) = -exp( eta(2,:)*bd2)
    A(6,5,1:np) = -exp(-eta(3,:)*bd2)*gamma(2,:)
    A(6,6,1:np) =  exp( eta(3,:)*bd2)*gamma(2,:)

    open(unit=44,file='debug_matrix.out',action='write',status='replace')
    do j=1,np
       if(j/=1) write(44,*) '  '
       write(44,'(A,I3,A,ES12.6,A,ES12.6,A)') 'p(',j,') = (',real(p(j)),',',aimag(p(j)),')'
       write(44,*) 'A:'
       do k=1,6
          write(44,565) A(k,:,j)
       end do
       write(44,*) 'b:'
       write(44,565) b(:,j)

       call zgesvx('EQ','NT', 6, 1, A(:,:,j), 6, AF(:,:,j), 6, &
            & ipiv(1:6), eqed, lpR, lpC, b(:,j), 6, nc(:,j), 6, rcond(j), ferr(j), &
            & berr(j),work,rwork,ierr)
       if(ierr /= 0) then
          write(44,*) 'ZGESV error code',ierr
       end if
       write(44,*) 'x:'
       write(44,565) nc(:,j)
       write(44,*) 'condition estimate of A:',rcond(j)
       write(44,*) 'forward error estimate of x:', ferr(j)  !! seems crazy
       write(44,*) 'backward error estimate of x:', berr(j)
    end do

565 format(6('(',ES11.3E3,',',ES11.3E3,')',1X))
566 format(6ES11.3E3)

    close(44)

    if(any(nc/=nc)) print *, '** NaN in matrix solution results **'

  end function numerical_matrix_twoaquifer
#endif
end module two_aquifer

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

module one_aquifer
  implicit none

  private
  public :: oneAquifer_leak

contains
  function oneAquifer_leak(dum) result(f)
    use shared_data, only :  lapalpha,laptol,tsval,lapM,zd,rd,piezometer, &
         &layer,obsGQy,obsGQw,obsGQz,alphaDz,alphaDr,alphaDy,Kdz,bd
    use constants, only : EP,DP,PI,RONE,RTWO,RZERO,PIEP,ETWO,EONE,EHALF
    use invlap, only : dehoog_invlap_ep    ! invlap.f90
    use complex_hyp_trig, only : cesinh, cecosh, cecoth,cesinhcosh
    use invlap_fcns

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none
    real(DP), intent(in) :: dum  ! scalar integration variable (Hankel parameter)
    real(DP) :: f, tee
    integer :: np, ord, k, np2, M2
    complex(EP), dimension(1:2*lapM+1) :: p, eta1,eta2,xi,gamma,delta,theis,fpep
    complex(EP), dimension(1:2*lapM+1) :: cosh1,sinh1,coth1
    complex(EP), dimension(4,1:2*lapM+1) :: c

    np = 2*lapM+1
    ord = size(obsGQy)
    call invlap_setup(p(1:np),tee)

    eta1(1:np) = sqrt((p(:) + real(dum,EP)**2)/real(alphaDz(1),EP))
    eta2(1:np) = sqrt((p(:) + real(alphaDr(2),EP)*real(dum,EP)**2)/real(alphaDz(2),EP))

    xi(1:np) = eta1(:)*real(alphaDy,EP)/p(:)
    !! the definition of gamma is different in old coordinates
    gamma(1:np) = eta1/(eta2*real(Kdz(1),EP))   
    theis(1:np) = ETWO/(p(:)*(p(:) + real(dum,EP)**2))

    coth1(1:np) = cecoth(eta2(:)*real(bd,EP))
    call cesinhcosh(eta1(:),sinh1,cosh1)

    !! this is the "new" delta, without the factor of 2
    delta(1:np) = (xi(:)*gamma(:)*coth1(:) + EONE)*sinh1(:) &
         & + (gamma(:)*coth1(:) + xi(:))*cosh1(:)

    ! observation is in unconfined pumped aquifer (factor of 2 is here now)
    if(layer(1)) then
       c(1,1:np) = -theis(:)/(ETWO*delta(:))* &
            & ((xi(:) - EONE)*exp(-eta1(:)) + gamma(:)*coth1(:) + EONE)
       c(2,1:np) = -theis(:)/(ETWO*delta(:))* &
            & ((xi(:) + EONE)*exp(eta1(:)) + gamma(:)*coth1(:) - EONE)

       if(piezometer) then
          fpep(1:np) = theis(:) + c(1,:)*exp(eta1(:)*real(zd,EP)) + &
               & c(2,:)*exp(-eta1(:)*real(zd,EP))
       else
          forall (k=1:ord)
             obsGQz(k,1:np) = theis(:) + c(1,:)*exp(eta1(:)*obsGQy(k)) + &
                  & c(2,:)*exp(-eta1(:)*obsGQy(k))
          end forall
          fpep(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
       end if

       ! observation is in confined unpumped aquitard
    elseif(layer(2)) then
       c(3,:) = 0.5_EP*theis(:)*(EONE - ((xi(:) - gamma(:))*cosh1(:) + &
            & (EONE - xi(:)*gamma(:))*sinh1(:) + gamma(:)*(coth1(:) + EONE))/delta(:))
       c(4,:) = 0.5_EP*theis(:)*(EONE - ((xi(:) + gamma(:))*cosh1(:) + &
            & (EONE + xi(:)*gamma(:))*sinh1(:) + gamma(:)*(coth1(:) - EONE))/delta(:))

       if(piezometer) then
          fpep(1:np) = c(3,:)*exp(eta2(:)*real(zd,EP)) + c(4,:)*exp(-eta2(:)*real(zd,EP))
       else
          forall (k=1:ord)
             obsGQz(k,1:np) = c(3,:)*exp(eta2(:)*obsGQy(k)) + c(4,:)*exp(-eta2(:)*obsGQy(k))
          end forall
          fpep(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
       end if

    else
       stop 'no layer 3 in one aquifer model'
    end if

    call check_nan(fpep,np2,M2)
    f = dum*dbesj0(dum*rd)* &
         & real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,fpep(1:np2),M2),DP)

  end function oneAquifer_leak
end module one_aquifer


!! ######################################################################
!! ######################################################################
module classical_mod
  implicit none

  private
  public :: classical_leak

contains

  function classical_leak(dum) result(f)
    use constants, only : EP,DP,RTWO,PIEP,EONE,ETWO 
    use shared_data, only : lapalpha,laptol,lapM,tsval,rd,zd,kappa,alphaDy,lambda,&
         & piezometer,obsGQy,obsGQz,obsGQw
    use invlap, only : dehoog_invlap_ep
    use complex_hyp_trig, only : cesinhcosh
    use invlap_fcns

#ifdef INTEL
    use ifport, only : dbesj0
#endif

    implicit none
    real(DP), intent(in) :: dum  ! scalar integration variable (Hankel parameter)
    real(DP) :: f, tee
    integer :: np, ord, k, np2, M2
    complex(EP), dimension(1:2*lapM+1) :: p,sinh1,cosh1,sinh2,cosh2
    complex(EP), dimension(1:2*lapM+1) :: theis, epfp, eta1,xi,omega,delta1

    np = 2*lapM+1
    ord = size(obsGQy)
    call invlap_setup(p(1:np),tee)

    eta1(1:np) = sqrt((p(:) + real(dum,EP)**2)/real(kappa,EP))
    xi(1:np) = eta1(:)*real(alphaDy,EP)/p(:)
    omega(1:np) = eta1(:)/real(lambda,DP)
    theis(1:np) = ETWO/(p(:)*(p(:) + real(dum,EP)**2))

    call cesinhcosh(eta1(:),sinh1,cosh1)
    delta1(1:np) = (omega(:) + xi(:))*cosh1(:) + (omega(:)*xi(:) + EONE)*sinh1(:)

    if(piezometer) then
       call cesinhcosh(eta1(:)*(1.0_EP - real(zd,EP)),sinh1,cosh1)
       call cesinhcosh(eta1(:)*real(zd,EP),sinh2,cosh2)

       epfp(1:np) = theis(:)*(EONE - &
            & (xi(:)*cosh1(:) + sinh1(:) + omega(:)*cosh2(:) + sinh2(:))/delta1(:))
    else
       do k=1,ord
          call cesinhcosh(eta1(:)*(1.0_EP - obsGQy(k)),sinh1,cosh1)
          call cesinhcosh(eta1(:)*obsGQy(k),sinh2,cosh2)

          obsGQz(k,1:np) = theis(:)*(EONE - &
               & (xi(:)*cosh1(:) + sinh1(:) + omega(:)*cosh2(:) + sinh2(:))/delta1(:))
       end do
       epfp(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
    end if

    call check_nan(epfp,np2,M2)   
    f = dum*dbesj0(dum*rd)*real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,epfp(1:np2),M2),DP)

  end function classical_leak
end module classical_mod

!! ######################################################################
!! ######################################################################
module neuman72_mod
  implicit none

  private
  public :: neuman72_noleak

contains

  function neuman72_noleak(dum) result(f)
    use constants, only : EP,DP,RTWO,ETWO,PIEP,EONE
    use shared_data, only : lapalpha,laptol,lapM,tsval,rd,zd,kappa,alphaDy,piezometer,&
         & obsGQy,obsGQw,obsGQz
    use invlap, only : dehoog_invlap_ep
    use complex_hyp_trig, only : cesinhcosh, cecosh
    use invlap_fcns

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none
    real(DP), intent(in) :: dum  ! scalar integration variable (Hankel parameter)
    real(DP) :: f, tee
    integer :: np, k, ord, np2, M2
    complex(EP), dimension(1:2*lapM+1) :: p,sinh1,cosh1
    complex(EP), dimension(1:2*lapM+1) :: theis,epfp,eta1,xi

    np = 2*lapM+1
    call invlap_setup(p(1:np),tee)

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

    if(piezometer) then 
       call cesinhcosh(eta1(:),sinh1,cosh1)
       ! one value of zd
       epfp(1:np) = theis(:)*(EONE - cecosh(eta1(:)*real(zd,EP))/ &
            & (xi(:)*sinh1(:) + cosh1(:)))
    else
       ! integrate over zd for each value of p using Gauss Quadrature
       ord = size(obsGQy)
       call cesinhcosh(eta1(:),sinh1,cosh1)
       do k=1,ord
          obsGQz(k,1:np) = theis(:)*(EONE - cecosh(eta1(:)*obsGqy(k))/ &
               & (xi(:)*sinh1(:) + cosh1(:)))
       end do
       ! (ztop-zbot) from Jacobian, 1/(ztop-zbot) from length => cancel 
       epfp(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
    end if

    call check_nan(epfp,np2,M2)
    f = dum*dbesj0(dum*rd)*real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,epfp(1:np2),M2),DP)

  end function neuman72_noleak
end module neuman72_mod


!! ######################################################################
!! ######################################################################
module neuman74_mod
  implicit none

  private
  public :: neuman74_noleak

contains

  function neuman74_noleak(dum) result(f)
    use constants, only : DP,EP,ETWO,EONE,EZERO
    use shared_data, only : lapalpha,laptol,lapM,tsval,rd,zd,ld,dd,alphaDz, &
         & alphaDy, obsGQy,obsGQw,obsGQz, piezometer
    use complex_hyp_trig, only : cecosh,cesinh,cesinhcosh
    use confined_partialpen, only : partialpen_laphank
    use invlap, only : dehoog_invlap_ep
    use invlap_fcns

#ifdef INTEL
    use ifport, only : dbesj0 
#endif

    implicit none

    real(DP), intent(in) :: dum
    real(EP) :: epld,epdd,epzd
    real(DP) :: f, tee
    complex(EP), dimension(0:1,1:2*lapM+1) :: ud
    complex(EP), dimension(1:2*lapM+1) :: fp, sinh1,cosh1, xi, eta1, p
    integer :: np, k, ord, np2, M2

    np = 2*lapM+1; epdd = real(dd,EP); epld = real(ld,EP); epzd = real(zd,EP)
    call invlap_setup(p(1:np),tee)
    eta1(1:np) = sqrt((p(:) + real(dum**2,EP))/real(alphaDz(1),EP))
    xi(1:np) = eta1(:)*real(alphaDy,EP)/p(:)
    call cesinhcosh(eta1(1:np),sinh1(1:np),cosh1(1:np))

    !! ud at watertable simplified for zd=0
    ud(0,1:np) = ETWO*(cesinh(eta1(:)*(EONE - epdd)) - &
         & cesinh(eta1(:)*(EONE - epld)))/&
         & (sinh1(:)*p(:)*eta1(:)**2*alphaDz(1)*(epld-epdd))

    if(piezometer) then 
       ud(1,1:np) = partialpen_laphank(eta1(1:np),epzd,p(1:np))

       ! one value of zd
       fp(1:np) = ud(1,:) - ud(0,:)*cecosh(eta1(:)*(EONE + epzd))/ &
            & (xi*sinh1(:) + cosh1(:))
    else
       ! integrate over zd for each value of p using Gauss Quadrature
       ord = size(obsGQy)
       do k=1,ord
          ud(1,1:np) = partialpen_laphank(eta1(1:np),obsGQy(k),p(1:np))
          obsGQz(k,1:np) = ud(1,:) - ud(0,:)*cecosh(eta1(:)*(EONE + obsGQy(k)))/ &
            & (xi*sinh1(:) + cosh1(:))
       end do
       ! (ztop-zbot) from Jacobian, 1/(ztop-zbot) from length => cancel 
       fp(1:np) = 0.5_EP*sum(spread(obsGQw(1:ord),dim=2,ncopies=np)*obsGQz(1:ord,1:np),dim=1)
    end if

    call check_nan(fp,np2,M2)
    f = dum*dbesj0(dum*rd)*real(deHoog_invlap_ep(lapalpha,laptol,tsval,tee,fp(1:np2),M2),DP)

  end function neuman74_noleak

end module neuman74_mod