! 
! 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
!
	
Program leakyUnconfinedDriver_version2
  ! $Id: driver2.f90,v 1.7 2007/04/03 00:21:37 kris Exp kris $

! needed for intel compiler
#ifdef INTEL
  use ifport, only : dbesj0, dbesj1     
#endif

  ! the finite portion of the inverse Hankel transform integral is evaluated
  ! using A. Miller's f90 implementation of the TOMS Algorithm 691 improvement
  ! on the dqags routine from the quadpack quadrature package.
  ! qxgs code downloaded from http://users.bigpond.net.au/amiller/NSWC/qxgs.f90  
  use adapt_quad, only : qxgs          ! qxgs.f90

  ! Pass data to routines using modules
  use shared_data         ! utility.f90

  ! constants and coefficients
  use constants, only : DP, SMALL, RONE, RZERO, PI
 
  ! all Laplace-Hankel space solutions are in one file (finite_laplace_hankel_solution.f90)
  ! each is its own module containing at least one function
  use neuman72_mod
  use neuman74_mod
  use confined_partialpen, only : partialpen_noleak
  use one_aquifer
  use two_aquifer
  use classical_mod

  use gauss_points

  ! approximation to exponential integral (for early and late theis curves)
  ! this is ~ single precision only (ok for plotting).
  use expint_approx

  implicit none

  real(DP), allocatable :: j0zero(:)
  real(DP) :: x, dx

  character(4) :: chint
  character(36) :: outfile, infile
  integer :: i, numt, terms, minlogsplit, maxlogsplit, splitrange, zerorange, obsintord
  logical :: oldcoord, lneuman72, lneuman74, lhantush64, lclassic, lambdacalc, l12Layer, lslimy
  real(DP), allocatable :: GLx(:),GLw(:)
  

  ! vectors of results and intermediate steps
  integer, allocatable :: splitv(:)
  real(DP), allocatable  :: ts(:), td(:), theisSs(:), theisSy(:)
  real(DP), allocatable :: finOneAqif(:), infOneAqif(:), oneAqif(:)
  real(DP), allocatable :: finTwoAqif(:), infTwoAqif(:), twoAqif(:)
  real(DP), allocatable :: finNeuman72(:), infNeuman72(:), neuman72(:)
  real(DP), allocatable :: finNeuman74(:), infNeuman74(:), neuman74(:)
  real(DP), allocatable :: finHantush64(:), infHantush64(:), hantush64(:)
  real(DP), allocatable :: finClassic(:), infClassic(:), classic(:)

  ! qxgs integration routine parameters
  real(DP) :: abserr, epsabs, epsrel1
  integer :: maxsub, ier, last, err, jj

  integer :: GLorder, j0split(1:2), naccel

  ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ! read input parameters from file some minimal error checking
  infile = 'input.dat'
  open(unit=19, file=infile, status='old', action='read', iostat=ier)
  if(ier /= 0) stop 'error opening input file for reading'

  read(19,*) oldcoord, l12Layer, piezometer, lneuman72, lneuman74, &
       & lhantush64, lclassic, lambdacalc, lslimy
  read(19,*) dd, ld
  read(19,*) b(1:MAXNLAYERS), lambda
  do i=1,MAXNLAYERS
     read(19,*) Kr(i), Kz(i)
  end do
  read(19,*) Ss(1:MAXNLAYERS)
  read(19,*) Sy
  read(19,*) numt
  read(19,*) zd, rd, ztop, zbot, obsintord

  if(any((/b(:),Kr(:),Kz(:),Ss(:),Sy,lambda/) < SMALL)) &
       & stop 'input error: zero or negative aquifer parameters'
  if(any((/Ss(:),Sy/) >= RONE)) stop 'input error: storage >= unity'
  
  if(.not. oldcoord) then
     bd2 = sum(b(1:2))/b(1)
     bd3 = sum(b(1:3))/b(1)
#ifdef DEBUG
     print *, 'bD2:',bd2
     print *, 'bD3:',bd3
#endif
  end if

  ! determine which layer the piezometer or wellscreen is in
  ! and perform some error checking on coordinates
  !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if(piezometer) then  ! one observation at zd (layer determined by zd)
     if (oldcoord) then  ! old coordinate system (+b1 > 0 > -b2)
        layer(3) = .false.
        if(zd > b(1) .or. zd < -b(2)/b(1)) then
           stop 'input error: z coordinate out of range (old system)'
        elseif(zd > -SMALL) then  ! unconfined aquifer
           layer(1) = .true.
           layer(2) = .false.
           print *, 'piezometer in aquifer'
        else                      ! aquitard (no flow bottom)
           layer(1) = .false.
           layer(2) = .true.
           print *, 'piezometer in aquitard'
        end if
     else ! new coordinate system  (0 > -b1 > -(b1+b2) > -(b1+b2+b3))
        if(zd > RZERO .or. zd < -bd3) then
           stop  'input error: z coordinate out of range (new system)'
        elseif(zd >= -RONE) then ! upper unconfined aquifer
           layer(1) = .true.
           layer(2) = .false.
           layer(3) = .false.
           print *, 'piezometer in upper aquifer'
        elseif(zd >= -bd2) then ! middle aquitard
           layer(1) = .false.
           layer(2) = .true.
           layer(3) = .false.
           print *, 'piezometer in aquitard'
        else                     ! lower confined aquifer
           layer(1) = .false.
           layer(2) = .false.
           layer(3) = .true.
           print *, 'piezometer in lower aquifer'
        end if
     end if
  else ! screened observation well (ztop & zbot used, zd not used)
     if(ztop <= zbot) stop 'z_top must be above z_bot'
     if(abs(ztop-zbot) < SMALL) stop 'z_top must be different from z_bot'
     if(oldcoord) then
        layer(3) = .false.
        if(ztop > b(1) .or. zbot < -b(2)/b(1)) stop 'z_bot or z_top out of range (old system)'
        ! both ends of observation wellscreen in aquifer
        if(ztop >= RZERO .and. zbot >= RZERO) then
           layer(1) = .true.
           layer(2) = .false.
           print *, 'screened obs well in aquifer'
        ! both ends of obervation wellscreen in aquitard
        elseif(ztop < RZERO .and. zbot < RZERO) then
           layer(1) = .false.
           layer(2) = .true.
           print *, 'screened obs well in aquitard'
        else
           stop 'ztop and zbot must both be in same unit (old coord system)'
        end if
     else
        if(ztop > RZERO .or. zbot < -bd3) stop 'z_bot or z_top out of range (new system)'
        !  wellscreen in upper unconfined aquifer
        if(ztop >= -RONE .and. zbot > -RONE) then
           layer(1) = .true.
           layer(2) = .false.
           layer(3) = .false.
           print *, 'screened obs well in upper aquifer'
        ! wellscreen in aquitard
        elseif(ztop >= -bd2 .and. zbot >= -bd2) then
           layer(1) = .false.
           layer(2) = .true.     
           layer(3) = .false.
           print *, 'screened obs well in aquitard'
        ! wellscreen in lower confined aquifer
        elseif(ztop >= -bd3 .and. zbot >= -bd3) then
           layer(1) = .false.
           layer(2) = .false.
           layer(3) = .true.
           print *, 'screened obs well in lower aquifer'
        else
           stop 'ztop and zbot must both be in same unit (new coord system)'
        end if
     end if ! oldcoord     
  end if ! piezometer

  if(oldcoord) then
     if(lneuman74 .or. lhantush64) then
        print *, '*** Neuman74 and Hantush64 not valid in old coord system **'
        lneuman74 = .false.
        lhantush64 = .false.
     end if
     if(layer(2) .and. (lclassic .or. lneuman72)) then
        print *, '*** classical leaky  and Neuman72 only valid in aquifer ***'
        lclassic = .false.
        lneuman72 = .false.
     end if
  else
     if(lneuman72 .or. lclassic) then
        print *,  '*** Neuman72 and Classical leaky invalid in new coord system **'
        lneuman72 = .false.
        lclassic = .false.
     end if
     if((layer(2) .or. layer(3)) .and. (lhantush64 .or. lneuman74)) then
        print *, '*** Hantush64 and Neuman74 only valid in upper aquifer ***'
        lhantush64 = .false.
        lneuman74 = .false.
     end if
     
     if(dd < RZERO .or. ld > RONE) stop 'dd or ld out of range (0-1)'
     if(dd > ld) stop 'ld must be > dd'
     if(abs(dd-ld) <= SMALL) stop 'dd and ld must be different'
  end if

  if(rd < SMALL) stop 'input error: rd must be > 0'
  read(19,*) lapalpha, laptol, lapM
  read(19,*) epsabs,epsrel1, maxsub, j0split(1:2), naccel, GLorder

  ! check if tolerances are too small (less than precision near 1)
  if(laptol < epsilon(1.0_DP)) then
     laptol = 10.0_DP*epsilon(1.0_DP)
     print *, 'tolerance for deHoog algorithm changed to: ',laptol
  end if
  if(epsabs < epsilon(1.0_DP)) then
     epsabs = 10.0_DP*epsilon(1.0_DP)
     print *, 'absolute tolerance for qxgs algorithm changed to: ', epsabs
  end if
  if(epsrel1 < epsilon(1.0_DP)) then
     epsrel1 = 10.0_DP*epsilon(1.0_DP)
     print *, 'relative tolerance for qxgs algorithm changed to: ', epsrel1
  end if

  if(any((/j0split(:),naccel,maxsub/) < 1)) &
       & stop 'max/min split, # accelerated terms and max # subroutine calls must be >= 1'
  ! solution vectors
  if(oldcoord) then
     allocate(oneAqif(numt), finOneAqif(numt), infOneAqif(numt))
     allocate(neuman72(numt), classic(numt))
     allocate(theisSs(numt), theisSy(numt))
     if(lneuman72) then
        allocate(finNeuman72(numt), infNeuman72(numt))
     end if
     if(lclassic) then
        allocate(finClassic(numt), infClassic(numt))
     end if
  else
     allocate(twoAqif(numt), finTwoAqif(numt), infTwoAqif(numt))
     allocate(neuman74(numt), hantush64(numt))
     if(lneuman74) then
        allocate(finNeuman74(numt), infNeuman74(numt))
     end if
     if(lhantush64) then
        allocate(finHantush64(numt), infHantush64(numt))
     end if
  end if
  
  allocate(ts(numt), td(numt), splitv(numt))
  terms = maxval(j0split(:))+naccel+1; allocate(j0zero(terms))
  read(19,*) ts(1:numt)
  read (19,*) outfile
  if(any(ts <= SMALL)) stop 'all times must be > 0'

  ! compute secondary aquifer properties (dimensionless)
  if(lambdacalc) then
     lambda = Kz(2)*b(1)/(b(2)*Kz(1))
  end if
  
  alphar(:) = Kr(:)/Ss(:)
  alphaz(:) = Kz(:)/Ss(:)
  Kdz(2:MAXNLAYERS) = Kz(2:MAXNLAYERS)/Kz(1:MAXNLAYERS-1)
  alphaDr(:) = alphar(:)/alphar(1)
  alphaDz(:) = alphaz(:)/alphar(1)
  bd = b(2)/b(1)    ! only used in old coord system calcs
  td(1:numt) = alphar(1)*ts(1:numt)/b(1)**2
  kappa = Kz(1)/Kr(1)  !! kappa == alphaDz(1), but kappa retained for old system
  alphaDy = alphaDz(1)*b(1)*Ss(1)/Sy

#ifdef DEBUG
  print *, 'lambda:',lambda
  print *, 'alphar:',alphar
  print *, 'alphaz:',alphaz
  print *, 'KDz:',Kdz
  print *, 'alphaDr:',alphaDr
  print *, 'alphaDz:',alphaDz
  print *, 'bd:',bd
  print *, 'kappa:',kappa
  print *, 'alphaDy:',alphaDy
#endif

  ! compute zeros of J0 bessel function
  !************************************
  ! asymptotic estimate of zeros - initial guess
  j0zero = (/((real(i-1,DP) + 0.75_DP)*PI,i=1,terms)/)
  do i=1,TERMS
     x = j0zero(i)
     NR: do jj = 1,100  
        ! Newton-Raphson
        dx = dbesj0(x)/dbesj1(x)
        x = x + dx
        if(abs(dx) < spacing(x)) then
           exit NR
        end if
     end do NR
     j0zero(i) = x
  end do

  ! split between finite/infinite part should be small for large time, large for small time
  do i=1,numt
     zerorange = j0split(2) - j0split(1)
     minlogsplit = floor(minval(log10(td)))   ! min log(td) -> maximum split
     maxlogsplit = ceiling(maxval(log10(td))) ! max log(td) -> minimum split
     splitrange = maxlogsplit - minlogsplit + 1 
     splitv(i) = j0split(1) + int(zerorange*((maxlogsplit - log10(td(i)))/splitrange))
  end do

#ifdef DEBUG
  print *, 'computed splits (integer zeros): ',splitv(:)
  print *, 'location of splits (zeros of J0):'
  do i= minval(splitv(:)),maxval(splitv(:))
     print *, i,j0zero(i)
  end do
#endif

  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ! finite portion of Hankel integral for each time level
  
  ! setup things for integrating over observation interval
  ! same between all functions

  if(.not. piezometer) then
     allocate(obsGQx(obsintord), obsGQw(obsintord), obsGQy(obsintord), &
          & obsGQz(obsintord,2*lapM+1))
     call gauss_quad_setup(obsintord,obsGQx,obsGQw)
     ! transform GQ abcissa to global z coordinate
     obsGQy(:) = 0.5_EP*(real(ztop - zbot,EP)*obsGQx(:) + real(ztop + zbot,EP))
     print *, 'Gauss Quadrature coordinates'
     print *, 'GQ space:',obsGQx
     print *, 'real space:',obsGQy
  end if

  write(*,*) 'integrate finite portion:'
  do i = 1, numt

     if (i == 1) then
        write (*,888) 'upper bound: ',j0zero(splitv(1))/rd
     elseif(splitv(i) /= splitv(i-1)) then
        write (*,888) 'upper bound: ',j0zero(splitv(i))/rd
     end if

     tsval = td(i)
     write (*,888) 'tD=', td(i)

     if(oldcoord) then

        ! previous two-layer system
        !===========================
        if(l12layer) then
           call qxgs(oneAquifer_leak,RZERO,j0zero(splitv(i))/rd,epsabs,epsrel1,&
                & finOneAqif(i),abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code two-layer: ', ier, abserr, last
           write(*,777) 'finOneAqif: ',finOneAqif(i)
        end if
        
        ! neuman 1972 
        ! ============
        if(lneuman72) then
           call qxgs(neuman72_noleak,RZERO,j0zero(splitv(i))/rd,epsabs,&
                & epsrel1,finNeuman72(i),abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code neuman72: ', ier, abserr, last
           write(*,777) 'finNeuman72: ',finNeuman72(i)
        end if

        ! classical leaky - unconfined 
        ! =============================
        if(lclassic) then
           call qxgs(classical_leak,RZERO,j0zero(splitv(i))/rd,epsabs,&
                & epsrel1,finClassic(i),abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code classical: ', ier, abserr, last
           write(*,777) 'finClassic: ',finClassic(i)
        end if

     else
        ! new three-layer system
        !===========================
        if(l12layer) then
           call qxgs(twoAquifer_leak,RZERO,j0zero(splitv(i))/rd,epsabs,epsrel1,finTwoAqif(i),&
                &abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code three-layer: ', ier, abserr, last
           write(*,777) 'finTwoAqif: ',finTwoAqif(i)
        end if
        
        ! neuman 1974 
        ! ============
        if(lneuman74) then
           call qxgs(neuman74_noleak,RZERO,j0zero(splitv(i))/rd,epsabs,&
                & epsrel1,finNeuman74(i),abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code neuman74: ', ier, abserr, last
           write(*,777) 'finNeuman74: ',finNeuman74(i)
        end if

        ! hantush 1964 
        ! =============
        if(lhantush64) then
           call qxgs(partialpen_noleak,RZERO,j0zero(splitv(i))/rd,epsabs,&
                & epsrel1,finHantush64(i),abserr, ier,maxsub, last)
           if (ier /= 0) write (*,*) 'QXGS: error code hantush64: ', ier, abserr, last
           write(*,777) 'finHantush64: ',finHantush64(i)
        end if
     end if

  end do

777 format(A15,ES14.6E3)
888 format(A,ES9.2)
  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ! "infinite" portion of Hankel  integral for each time level
  ! integrate between zeros of Bessel function, extrapolate 
  ! area from series of areas of alternating sign

  write (*,*) 'integrate infinite portion'
  allocate(Glx(GLorder-2),GLw(GLorder-2))
  call gauss_lobatto_setup(GLorder,GLx,GLw)  ! get weights and abcissa
  do i = 1, numt

     if (i == 1) then
        write(*,888) 'lower bound:', j0zero(splitv(1))/rd
     elseif(splitv(i) /= splitv(i-1)) then
        write(*,888) 'lower bound:', j0zero(splitv(i))/rd
     end if
     tsval = td(i)
     write (*,888) 'tD=', td(i)

     if(oldcoord) then
        ! two-layer system
        !===========================
        if(l12layer) then
           infOneAqif(i) = inf_integral(oneAquifer_leak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infOneAqif: ',infOneAqif(i)
        end if
        
        ! neuman 1972 
        !=============
        if(lneuman72) then           
           infNeuman72(i) = inf_integral(neuman72_noleak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infNeuman72: ',infNeuman72(i)
        end if

        ! "classical" leaky-unconfined 
        !============================
        if(lclassic) then
           infClassic(i) = inf_integral(classical_leak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infClassic: ',infClassic(i)
        end if

     else
        ! new three-layer system
        !===========================
        if(l12layer) then
           infTwoAqif(i) = inf_integral(twoAquifer_leak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infTwoAqif: ',infTwoAqif(i)
        end if
        
        ! neuman 1974 
        !=============
        if(lneuman74) then           
           infNeuman74(i) = inf_integral(neuman74_noleak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infNeuman74: ',infNeuman74(i)
        end if

        ! hantush 1964
        !==============
        if(lhantush64) then
           infHantush64(i) = inf_integral(partialpen_noleak,splitv,naccel,glorder,j0zero,GLx,GLw)
           write(*,777) 'infHantush64: ',infHantush64(i)
        end if
     end if
  end do
  

  !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ! Theis and Neuman solutions (for comparison)
  if(oldcoord) then
     do i = 1, numt
        theisSs(i) = expint(rD*rD/(4.0_DP*tD(i)))             ! early Theis (Ss)
        theisSy(i) = expint(rD*rD*Sy*b(1)/(4.0_DP*Kz(1)*ts(i)))  ! late Theis (Sy)
     end do
  end if
  
  ! add parts together
  if(oldcoord) then
     
     !! silently throw away infinite portion if it didn't converge (lslimy == true)
     if(lslimy) then
        where(abs(infOneAqif(:)) < huge(1.0_DP))
           oneAqif(1:numt) = finOneAqif(:) + infOneAqif(:)
        elsewhere
           oneAqif(1:numt) = finOneAqif(:)
        end where
     else
        oneAqif(1:numt) = finOneAqif(:) + infOneAqif(:)
     end if
     
     if(lneuman72) then
        if(lslimy) then
           where(abs(infNeuman72) < huge(1.0_DP))
              neuman72(1:numt) = finNeuman72(:) + infNeuman72(:)
           elsewhere
              neuman72(1:numt) = finNeuman72(:)
           end where
        else
           neuman72(1:numt) = finNeuman72(:) + infNeuman72(:)
        end if
     else
        neuman72(:) = RZERO
     end if
     
     if(lclassic) then
        if(lslimy) then
           where(abs(infClassic) < huge(1.0_DP))
              classic(1:numt) = finClassic(:) + infClassic(:)
           elsewhere
              classic(1:numt) = finClassic(:)
           end where
        else
           classic(1:numt) = finClassic(:) + infClassic(:)
        end if
     else
        classic(:) = RZERO
     end if
     
  else !! new coordinate system
     if(lslimy) then
        where(abs(infTwoAqif) < huge(1.0_DP))
           twoAqif(1:numt) = finTwoAqif(:) + infTwoAqif(:)
        elsewhere
           twoAqif(1:numt) = finTwoAqif(:)
        end where
     else
        twoAqif(1:numt) = finTwoAqif(:) + infTwoAqif(:)
     end if

     if(lneuman74)  then
        if(lslimy) then
           where(abs(infNeuman74) < huge(1.0_DP))
              neuman74(1:numt) = finNeuman74(:) + infNeuman74(:)
           elsewhere
              neuman74(1:numt) = finNeuman74(:)
           end where
        else
           neuman74(1:numt) = finNeuman74(:) + infNeuman74(:)
        end if
     else
        neuman74(:) = RZERO
     end if
     
     if(lhantush64) then
        if(lslimy) then
           where(abs(infHantush64) < huge(1.0_DP))
              hantush64(1:numt) = finHantush64(:) + infHantush64(:)
           elsewhere
              hantush64(1:numt) = finHantush64(:)
           end where
        else
           hantush64(1:numt) = finHantush64(:) + infHantush64(:)
        end if
     else
        hantush64(:) = RZERO
     end if
     
  end if

  ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  ! write results to file 

  ! open output file or die
  open (unit=20, file=outfile, status='replace', action='write', iostat=err)
  if (err /= 0) stop 'cannot open output file for writing'

  ! echo input parameters
  write (20,'(6(A,L2))') '# old coordinates:',oldcoord,' pt soln:', piezometer, &
       & ' Neuman72 soln:',lneuman72, ' classical soln:', lclassic, ' Neuman74 soln:',&
       & lneuman74, ' Hantush64 soln:',lhantush64 
  write (20,'(A,3(1X,ES10.3))') '# b:',b(1:3)
  write (20,'(A,3(1X,ES10.3))') '# Kr:',Kr(1:3)
  write (20,'(A,3(1X,ES10.3))') '# Kz:',Kz(1:3)
  write (20,'(A,3(ES10.3),A,ES10.3)') '# Ss:',Ss(1:3), ' Sy:',Sy
  if (piezometer) then
     write (20,'(2(A,ES10.3))') '# obs location zd:',zd,' rd:', rd
  else
     write (20,'(2(A,ES10.3),A,I3)') '# obs interval ztop:',ztop,'  zbot:',zbot,&
          & ' Gauss quad order:',obsintord
  end if
  
  write (20,'(2(A,ES10.3),A,I4)') '# L^{-1} parameters alpha:',lapalpha,' tol:',laptol,' M:',lapM
  write (20,'(2(A,ES10.3),A,I6)') '# H^{-1} qxgs params  epsabs:', epsabs, ' epsrel',epsrel1, &
       & ' max subs:',maxsub
  write (20,'(A,2I3,ES10.3,A,I4)') '# min/max J0(x) zero split:', j0split(1:2), &
       & j0zero(splitv(1)), ' # 0s to accelerate:', naccel
  write (20,'(A,I2)') '# order of Gauss-Lobatto quad:', GLorder 
  write(chint,'(I4.4)') size(ts,1)
  write (20,'(A,'//chint//'(1X,ES10.4))') '# original (dimensional) time:', ts
  write (20,'(A)') '#'
  
  if(oldcoord) then
     write (20,'(A)') '#         t_D                Theis (Ss)                Theis (Sy)      &
          &          s_{D} classic               one_layer'
     write (20,'(A)') '#----------------------------------------------------------------------&
          &--------------------------------------------------'
     do i = 1, numt
        write (20,'(5(1x,ES24.15E3))') td(i),theisSs(i),theisSy(i),classic(i),oneAqif(i)
     end do

     if(lneuman72) then
        write (20,*) ' '
        write (20,*) ' '
        write (20,'(A)') '#         t_D              Neuman72  '
        write (20,'(A)') '#------------------------------------'
        do i = 1, numt
           write (20,'(2(1x,ES24.15E3))') td(i),neuman72(i)
        end do
  end if
     
  else  ! new coord system
     write (20,'(A)') '#         t_D                two_layer'
     write (20,'(A)') '#--------------------------------------'
     do i = 1, numt
        write (20,'(2(1x,ES24.15E3))') td(i),twoAqif(i)
     end do

     if(lhantush64 .or. lneuman74) then
        write (20,*) ' '
        write (20,*) ' '
        write (20,'(A)') '#         t_D              Hantush64               Neuman74 '
        write (20,'(A)') '#-------------------------------------------------------------'
        do i = 1, numt
           write (20,'(4(1x,ES24.15E3))') td(i),hantush64(i),neuman74(i)
     end do
  end if
end if

contains

  !! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  
  !!   end of program flow
  !! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  

  !! ###################################################
  function inf_integral(integrand,split,num,ord,zeros,GLx,GLw) result(integral)
    use constants, only : DP, HALF
    use shared_data, only : rd

    interface 
       function integrand(a) result(H)
         use constants, only : DP
         real(DP), intent(in) :: a
         real(DP) :: H
       end function integrand
    end interface

    integer, dimension(:), intent(in) :: split
    integer, intent(in) :: num,ord
    real(DP), dimension(:), intent(in) :: zeros

    real(DP), dimension(ord-2) :: Gly,Glz
    real(DP), dimension(:), intent(in) :: GLx,GLw
    real(DP), dimension(num) :: area
    real(DP) :: lob,hib,width,integral
    integer :: j,k

    do j = split(i)+1, split(i)+num
       lob = zeros(j-1)/rd ! lower bound
       hib = zeros(j)/rd   ! upper bound
       width = hib - lob
       
       ! transform GL abcissa to global coordinates
       Gly(:) = HALF*(width*GLx(:) + (hib + lob))
       
       do k = 1,ord-2
          GLz(k) = integrand(Gly(k))
       end do
       area(j-split(i)) = HALF*width*sum(GLz(1:ord-2)*GLw(1:ord-2))
    end do

    integral = wynn_epsilon(area(1:num))

  end function inf_integral

  !! ###################################################
  subroutine gauss_lobatto_setup(ord,abcissa,weight)
    integer, intent(in) :: ord
    real(DP), intent(out), dimension(ord-2) :: weight, abcissa

    ! leave out the endpoints (abcissa = +1 & -1), since 
    ! they will be the zeros of the Bessel functions
    ! (therefore, e.g., 5th order integration only uses 3 points)
    ! copied from output of Matlab routine by Greg von Winckel, at
    ! http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4775

    select case(ord)
    case(5)
       abcissa(1:2) = (/-0.65465367070798_DP, 0.0_DP/)
       abcissa(3) = -abcissa(1)
       weight(1:2) =   (/0.54444444444444_DP, 0.71111111111111_DP/)
       weight(3) = weight(1)
    case(6)
       abcissa(1:2) = (/-0.76505532392946_DP, -0.28523151648065_DP/)
       abcissa(3:4) = -abcissa(2:1:-1)
       weight(1:2) = (/0.37847495629785_DP, 0.55485837703549_DP/)
       weight(3:4) = weight(2:1:-1)
    case(7)
       abcissa(1:3) = (/-0.83022389627857_DP, -0.46884879347071_DP, 0.0_DP/)
       abcissa(4:5) = -abcissa(2:1:-1)
       weight(1:3) = (/0.27682604736157_DP, 0.43174538120986_DP, &
            & 0.48761904761905_DP/)
       weight(4:5) = weight(2:1:-1)
    case(8)
       abcissa(1:3) = (/-0.87174014850961_DP, -0.59170018143314_DP, &
            & -0.20929921790248_DP/)
       abcissa(4:6) = -abcissa(3:1:-1)
       weight(1:3) = (/0.21070422714351_DP, 0.34112269248350_DP, &
            & 0.41245879465870_DP/)
       weight(4:6) = weight(3:1:-1)
    case(9)
       abcissa(1:4) = (/-0.89975799541146_DP, -0.67718627951074_DP, &
                      & -0.36311746382618_DP, 0.0_DP/)
       abcissa(5:7) = -abcissa(3:1:-1)
       weight(1:4) = (/0.16549536156081_DP, 0.27453871250016_DP, &
            & 0.34642851097305_DP, 0.37151927437642_DP/)
       weight(5:7) = weight(3:1:-1)
    case(10)
       abcissa(1:4) = (/-0.91953390816646_DP, -0.73877386510551_DP, &
            & -0.47792494981044_DP, -0.16527895766639_DP/)
       abcissa(5:8) = -abcissa(4:1:-1)
       weight(1:4) = (/0.13330599085107_DP, 0.22488934206313_DP, &
            & 0.29204268367968_DP, 0.32753976118390_DP/)
       weight(5:8) = weight(4:1:-1)
    case(11)
       abcissa(1:5) = (/-0.93400143040806_DP, -0.78448347366314_DP, &
            & -0.56523532699621_DP, -0.29575813558694_DP, 0.0_DP/)
       abcissa(6:9) = -abcissa(4:1:-1)
       weight(1:5) = (/0.10961227326699_DP, 0.18716988178031_DP, &
            & 0.24804810426403_DP, 0.28687912477901_DP, &
            & 0.30021759545569_DP/)
       weight(6:9) = weight(4:1:-1)
    case(12)
       abcissa(1:5) = (/-0.94489927222288_DP, -0.81927932164401_DP, &
            & -0.63287615303186_DP, -0.39953094096535_DP, &
            & -0.13655293285493_DP/)
       abcissa(6:10) = -abcissa(5:1:-1)
       weight(1:5) = (/0.09168451741320_DP, 0.15797470556437_DP, &
            & 0.21250841776102_DP, 0.25127560319920_DP, &
            & 0.27140524091070_DP/) 
       weight(6:10) = weight(5:1:-1)
    case default
       stop 'unsupported Gauss-Lobatto integration order selected'
    end select
  end subroutine gauss_lobatto_setup


  !! ###################################################
  !! wynn-epsilon acceleration of partial sums, given a series
  function wynn_epsilon(series) result(accsum)
    use constants, only : DP, RONE, RZERO

    integer, parameter :: MINTERMS = 4
    real(DP), dimension(1:), intent(in) :: series
    real(DP) :: denom, accsum
    integer :: np, i, j, m

    real(DP), dimension(1:size(series),-1:size(series)-1) :: ep

    np = size(series)

    ! first column is intiallized to zero
    ep(:,-1) = RZERO

    ! build up partial sums, but check for problems
    check: do i=1,np
       if((abs(series(i)) > huge(RZERO)).or.(series(i) /= series(i))) then
          ! +/- Infinity or NaN ? truncate series to avoid corrupting accelerated sum
          np = i-1
          if(np < MINTERMS) then
             write(*,*) 'not enough useful terms in Wynn-Epsilon series to accelerate'
             accsum = 1.0/0.0  ! make it clear answer is bogus
             goto 777
          else
             write(*,'(A,I3,A)') 'Wynn-Epsilon series&
                  &, truncated to ',np,' terms.'
             exit check
          end if
       else
          ! term is good, continue
          ep(i,0) = sum(series(1:i))
       end if
    end do check
    
    ! build up epsilon table (each column has one less entry)
    do j=0,np-2 
       do m = 1,np-(j+1)
          denom = ep(m+1,j) - ep(m,j)
          if(abs(denom) > tiny(1.0_DP)) then ! check for div by zero
             ep(m,j+1) = ep(m+1,j-1) + RONE/denom
          else
             accsum = ep(m+1,j)
             write(*,'(A,I2,1X,I2)') 'cancellation in epsilon table ',m,j
             goto 777
          end if
       end do
    end do

    ! if made all the way through table use "corner value" of triangle as answer
    if(mod(np,2) == 0) then
       accsum = ep(2,np-2)  ! even number of terms - corner is acclerated value
    else
       accsum = ep(2,np-3)  ! odd numbers, use one column in from corner
    end if
777 continue

  end function wynn_epsilon

end program leakyUnconfinedDriver_version2
