! ! 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 ! 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