! ! Copyright (c) 2008 Kristopher L. Kuhlman (kuhlman at hwr dot arizona dot edu) ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the Software is ! furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN ! THE SOFTWARE. ! ! Malama, B., A. Revil, and K.L. Kuhlman, 2009. A semi-analytical solution ! for transient streaming potentials associated with confined aquifer ! pumping tests, Geophysical Journal International, 176(3), 1007–1016. ! http://dx.doi.org/10.1111/j.1365-246X.2008.04014.x ! ! Malama, B., K.L. Kuhlman, and A. Revil, 2009. Theory of transient ! streaming potentials associated with axial-symmetric flow in unconfined ! aquifers, Geophysical Journal International, 179(2), 990–1003. ! http://dx.doi.org/10.1111/j.1365-246X.2009.04336.x ! Program Driver ! $Id: driver4.f90,v 1.2 2009/03/22 20:06:54 kris Exp kris $ ! needed for intel compiler #ifdef INTEL use ifport, only : dbesj0, dbesj1 #endif ! Pass data to routines using modules use shared_data ! utility.f90 ! constants and coefficients use constants, only : DP, SMALL, RONE, RZERO, PI use invlap_fcns, only : stehfest_coeff !! three-layer solution use three_layer implicit none real(DP), allocatable :: j0zero(:) real(DP) :: x, dx character(23) :: fmt character(75) :: outfile, infile integer :: i, j, w, numt, terms, minlogsplit, maxlogsplit integer :: splitrange, zerorange real(DP), allocatable :: GLx(:),GLw(:) ! vectors of results and intermediate steps integer, allocatable :: splitv(:) real(DP), allocatable :: ts(:), td(:), dt(:) , fa(:) real(DP), allocatable :: finOneAqif(:,:), infOneAqif(:,:), totOneAqif(:,:) real(DP), allocatable :: oneAqifD(:), deriv(:) ! coordinates of pumping, injecting, image, and observation wells ! general reflection matrix ! distances from observation point to real & image wells real(DP), allocatable :: coord(:,:), r(:) real(DP), dimension(2) :: obs real(DP) :: depth, distb integer :: npts,nwells,BC, nn logical :: Lbdry !! k parameter in tanh-sinh quadrature (order is 2^k - 1) integer :: tsk, tsN, nst real(DP), allocatable :: tsw(:), tsa(:), sgn(:), tshh(:), tmp(:) integer, allocatable :: tskk(:), tsNN(:), ii(:) integer :: GLorder, j0split(1:2), naccel, err, m ! realted to dimensional solution Logical :: quiet, dipole real(DP) :: HC, PhiC, Qrate, PolErr ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! read input parameters from file some minimal error checking infile = 'input.dat' open(unit=19, file=infile, status='old', action='read') ! suppress output to screen, T/F dipole test, integer BC type, ! perpendicular distance to boundary from pumping well read(19,*) quiet, dipole, BC, distb ! dipole test has 2 pumping wells if (dipole) then nwells = 2 else nwells = 1 end if ! not actually used if not dipole test (pumping well at origin automatically) allocate(coord(2,nwells)) if (BC > 2 .or. BC < 0) then print *, 'BC = 0 for no boundary, 1 constant head, 2 for no-flow' stop end if coord(1:2,1) = (/ 0.0,0.0 /) ! pumping well coordinates if (nwells == 1) then read(19,*) read(19,*) else read(19,*) coord(1,2:nwells) ! read in x-coordinates of pumping wells read(19,*) coord(2,2:nwells) ! "" y- "" "" end if if (BC == 1 .or. BC == 2) then ! if constant head (1) or no-flow (2) BC, need image wells Lbdry = .true. allocate(r(2*nwells),sgn(2*nwells)) else Lbdry = .false. allocate(r(nwells),sgn(nwells)) end if read(19,*) Qrate ! injection is -Qrate, pumping is +Qrate read(19,*) b(1:3) ! layer thicknesses read(19,*) Kr, kappa, Ss, Sy read(19,*) sigma(1:3) read(19,*) gamell(1:3) !! rho*g*ell = gamma*ell = C (only used to re-dimensionalize) read(19,*) numt read(19,*) obs(1:2),depth ! x,y,depth of observation point if (.not. quiet) then write(*,'(2(A,L1),A,I1,A,ES10.3)') 'quiet?', quiet, ' dipole?',dipole, ' BC:',BC, ' dist Q -> bdry:',distb write(*,111) 'layer b:',b(1:3) write(*,'(A,4(1X,ES10.3))') 'Kr,kappa,Ss,Sy:', Kr, kappa, Ss, Sy write(*,111) 'sigma:',sigma(1:3) write(*,111) 'gamma*ell:',gamell(1:3) write(*,'(A,1X,I4)') '# times',numt write(*,222) '(x_Q,y_Q): (',coord(1,1),', ',coord(2,1), ')' if (dipole) then write(*,222) '(x_I,y_I: (',coord(1,2),', ',coord(2,2),')' end if write(*,222) '(x_obs,y_obs): (',obs(1),', ',obs(2),')' write(*,'(A,ES10.3)') 'obs depth:',depth end if 111 format(A,3(1X,ES10.3)) 222 format(2(A,ES10.3),A) ! need some error-checking for electrical parameters if(any((/b(:),Kr,kappa,Ss,Sy,sigma(:)/) < SMALL)) then print *, 'input error: zero or negative aquifer parameters' stop end if if(distb <= 0.0) then print *, 'distb should be positive' stop end if if(obs(1) < -distb) then print *, 'observation point must be on this side of boundary' stop end if if (dipole) then if (coord(1,2) < -distb) then print *, 'injection point must be on this side of boundary' stop end if end if ! dimensionless thicknesses bd1 = b(1)/b(2) bd3 = b(3)/b(2) sgn(:) = -999.99 ! Q -> obs r(1) = sqrt(obs(1)**2 + obs(2)**2) sgn(1) = +1.0 if (.not. quiet) write(*,*) '---------------------------' if (Lbdry) then if (dipole) then ! I -> obs r(2) = sqrt((obs(1)-coord(1,2))**2 + (obs(2)-coord(2,2))**2) sgn(2) = -1.0 ! I image -> obs r(4) = sqrt((obs(1)+2.0_DP*distb+coord(1,2))**2 + (obs(2)-coord(2,2))**2) if (BC == 2) then sgn(4) = -1.0 elseif (BC == 1) then sgn(4) = +1.0 end if ! Q image -> obs r(3) = sqrt((obs(1)+2.0_DP*distb)**2 + obs(2)**2) if (BC == 2) then sgn(3) = +1.0 elseif (BC == 1) then sgn(3) = -1.0 end if npts = 4 if (.not. quiet) then write(*,'(A)') ' '//' Q '//' Inj '//& & ' Q img '//' I Img ' end if else ! Q image -> obs r(2) = sqrt((obs(1)+2.0_DP*distb)**2 + obs(2)**2) if (BC == 2) then sgn(2) = +1.0 elseif (BC == 1) then sgn(2) = -1.0 end if npts = 2 if (.not. quiet) then write(*,'(A)') ' '//' Q '//' Q img ' end if end if else if (dipole) then ! I -> obs r(2) = sqrt((obs(1)-coord(1,2))**2 + (obs(2)-coord(2,2))**2) sgn(2) = -1.0 npts = 2 if (.not. quiet) then write(*,'(A)') ' '//' Q '//' I ' end if else npts = 1 if (.not. quiet) then write(*,'(A)') ' '//' Q ' end if end if end if ! non-dimensionalize all distances by aquifer thickness (second layer) r(:) = r(:)/b(2) ! zd is positive up from bottom of aquifer; depth down from surface zd = 1.0_DP + (b(1) - depth)/b(2) if (.not. quiet) then fmt = '(A, (ES10.3,1X)) ' write(fmt(4:4),'(I1.1)') npts write(*,fmt) 'r_D: ',r fmt(5:15) = '(1I10,001X)' write(*,fmt) 'sign of Q:',nint(sgn) write(*,*) '---------------------------' write(*,'(A,ES10.3)') 'z_D: ',zd end if ! perform some error checking on coordinates !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% layer = .false. if(zd > bd1 + 1.0 .or. zd < -bd3) then print *, 'input error: z coordinate out of range' stop elseif(zd >= 1.0) then ! upper unpumped vadose zone layer(1) = .true. if(.not. quiet) then print *, 'observation point in vadose zone above unconfined layer' end if elseif(zd > 0.0) then ! middle pumped unconfined aquifer layer(2) = .true. if(.not. quiet) then print *, 'observation point in pumped unconfined aquifer' end if else ! lower unpumped aquiclude layer(3) = .true. if(.not. quiet) then print *, 'observation point in lower unpumped aquiclude' end if end if read(19,*) stehfestN allocate(stehfestV(stehfestN)) stehfestV(:) = stehfest_coeff(stehfestN) read(19,*) tsk, j0split(1:2), naccel, GLorder, nst if(tsk - nst < 2) then print *, 'Tanh-Sinh k is too low (',tsk,') for given level& & of Richardson extrapolation (',nst,'). Increase tsk or decrease nst.' stop end if if(any((/j0split(:),naccel, tsk/) < 1)) then print *, 'max/min split, # accelerated terms and tanh-sinh k must be >= 1' stop end if ! solution vectors allocate(oneAqifD(numt), finOneAqif(numt,npts), infOneAqif(numt,npts)) allocate(deriv(numt), totOneAqif(numt,npts)) allocate(ts(numt), td(numt), splitv(numt), dt(numt)) !! ,ud(numt)) terms = maxval(j0split(:))+naccel+1 allocate(j0zero(terms)) ! time is now read as a column rather than a long row, ! to accomidate the fact that PEST has a maximum row length do i=1,numt read(19,*) ts(i) end do read (19,*) outfile if(any(ts <= SMALL)) then print *, 'all times must be > 0' stop end if alphaR = Kr/Ss Kz = kappa*Kr alphaY = b(2)*Kz/Sy alphaD = alphaY/alphaR sigmaD(1:3) = [ sigma(1)/sigma(2),1.0_DP,sigma(3)/sigma(2) ] td(1:numt) = alphaR*ts(1:numt)/b(2)**2 #ifdef DEBUG write(*,333) 'alpha_r:',alphaR write(*,333) 'K_z:', Kz write(*,333) 'alpha_y:',alphaY write(*,333) 'alpha_D:',alphaD write(*,'(A,3(1X,ES10.3))') 'sigma_D:',sigmaD write(*,333) 'b_{D,1}:',bd1 write(*,333) 'b_{D,3}:',bd3 333 format(A,ES10.3) #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 ! 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 = maxval(j0split(:)) - minval(j0split(:)) 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) = minval(j0split(:)) + & & 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 !!$ open(unit=222,file='debug.dbg',action='write',status='replace') do w = 1,npts rD = r(w) if (.not. quiet) then write(*,'(A,I3,A,ES10.3)') 'integrating obs -> well ',w,' rD:',rD end if !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! finite portion of Hankel integral for each time level if(.not. quiet) then write(*,*) 'integrate finite portion:' write(*,*) '==========================' end if if(.not. allocated(tsw)) then tsN = 2**tsk - 1 allocate(tsw(tsN),tsa(tsN),fa(tsN),ii(tsN)) allocate(tskk(nst),tsNN(nst),tshh(nst),tmp(nst)) end if do i = 1, numt if(.not. quiet) then 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 end if tsval = td(i) ! three-layer unconfined electrokinetic system !=========================== tskk(1:nst) = (/ (tsk-m,m=nst-1,0,-1)/) tsNN(1:nst) = 2**tskk - 1 tshh(1:nst) = 4.0_DP/(2**tskk) ! compute abcissa and for highest-order case (others are subsets) call tanh_sinh_setup(tskk(nst),j0zero(splitv(i))/rD,tsw(1:tsNN(nst)),tsa(1:tsNN(nst))) ! compute solution at densest set of abcissa do nn=1,tsNN(nst) fa(nn) = unconfined_aquifer_sp( tsa(nn) ) end do tmp(nst) = (j0zero(splitv(i))/rD)/2.0_DP*sum(tsw(:)*fa(:)) do j=nst-1,1,-1 ! only need to re-compute weights for each subsequent step call tanh_sinh_setup(tskk(j),j0zero(splitv(i))/rD,tsw(1:tsNN(j)),tsa(1:tsNN(j))) ! compute index vector, to slice up solution ! for nst'th turn count regular integers ! for (nst-1)'th turn count even integers ! for (nst-2)'th turn count every 4th integer, etc... ii(1:tsNN(j)) = (/( m* 2**(nst-j), m=1,tsNN(j) )/) ! estimate integral with subset of function evaluations and appropriate weights tmp(j) = (j0zero(splitv(i))/rD)/2.0_DP*sum(tsw(1:tsNN(j))*fa(ii(1:tsNN(j)))) end do if (nst > 1) then call polint(tshh(1:nst),tmp(1:nst),0.0_DP,finOneAqif(i,w),PolErr) else finOneAqif(i,w) = tmp(1) end if if(.not. quiet) then write(*,777) 'tD=', td(i),' ',finOneAqif(i,w) end if end do 777 format(A,ES8.2,A,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 if(.not. quiet) then write(*,*) write (*,*) 'integrate infinite portion' write(*,*) '=============================' end if if(.not. allocated(GLx)) then allocate(Glx(GLorder-2),GLw(GLorder-2)) call gauss_lobatto_setup(GLorder,GLx,GLw) ! get weights and abcissa end if do i = 1, numt if(.not. quiet) then 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 end if tsval = td(i) ! new three-layer unconfined system !=========================== infOneAqif(i,w) = inf_integral(unconfined_aquifer_sp,splitv,naccel,& &glorder,j0zero,GLx,GLw) if(.not. quiet) then write(*,777) 'tD=', td(i),' ',infOneAqif(i,w) end if end do end do ! combine solutions and re-dimensionalize HC = Qrate/(4.0_DP*PI*b(2)*Kr) PhiC = gamell(2)*HC/sigma(2) totOneAqif(1:numt,1:npts) = (finOneAqif(:,:) + infOneAqif(:,:))*PhiC oneAqifD = 0.0_DP ! just add all wells (real and image), sign of wells is handled above ! use npts, rather than nwells to handle both dipole and no dipole cases do j=1,npts oneAqifD(1:numt) = oneAqifD + sgn(j)*totOneAqif(:,j) end do ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! compute logarithmic deriviative ds/d(ln(t)) = t*(ds/dt) ! using finite-differences in time ! use backwards diff for all but first step deriv(2:) = ts(2:)*(oneAqifD(2:) - oneAqifD(1:))/(ts(2:) - ts(1:numt-1)) deriv(1) = ts(1)*(oneAqifD(1) - oneAqifD(2))/(ts(1) - ts(2)) ! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ! write results to file ! open output file or die open (unit=20, file=outfile, status='replace', action='write', iostat=err) if (err /= 0) then print *, 'cannot open output file for writing' stop end if ! echo input parameters write (20,'(1(A,L2))') '# suppress output?:',quiet write (20,'(A,1(ES10.3,1X))') '# pumping rate Q', Qrate write (20,'(A,3(1X,ES10.3))') '# b:',b(1:3) write (20,'(A,3(1X,ES10.3))') '# sigma:',sigma(1:3) write (20,'(A,3(1X,ES10.3))') '# gamma*ell:',gamell(1:3) write (20,'(2(A,ES10.3))') '# Ss:',Ss, ' Kr:',Kr write (20,'(2(A,ES10.3))') '# Sy:',Sy, ' Kz:',Kz write (20,'(3(A,ES10.3))') '# obs zd:',zd fmt = '(A,I2,A, (ES10.3,1X))' write(fmt(9:10),'(I2.2)') stehfestN write (20,fmt) '# L^{-1} parameter N:',stehfestN, ' V:', stehfestV write (20,'(A,I4)') '# tanh-sinh quadrature order:',tsN write (20,'(A,2I3,ES10.3,A,I4)') & &'# min/max J0(x) zero split:', minval(j0split(:)),maxval(j0split(:)), & & j0zero(splitv(1)), ' # 0s to accelerate:', naccel write (20,'(A,I2)') '# order of Gauss-Lobatto quad:', GLorder write (20,'(A)') '#' write (20,'(A)') '# t_D & & H^{-1}[ L^{-1}[ f_D ]] & & log-deriv ' write (20,'(A)') '#----------------------------& &------------------------------------------------' do i = 1, numt write (20,'(3(1x,ES24.15E3))') ts(i), oneAqifD(i), deriv(i) end do deallocate(stehfestV) 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 implicit none 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 #ifdef DEBUG if (any(GLz /= GLz)) then print *, 'NaN in GLz' elseif (any(abs(GLz) >= huge(1.0D1))) then print *, '+/- infinity in GLz' end if #endif 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 tanh_sinh_setup(k,s,w,a) use constants, only : PIOV2 implicit none integer, intent(in) :: k real(DP), intent(in) :: s real(DP), intent(out), dimension(2**k-1) :: w, a integer :: N,r,i real(DP) :: h real(DP), allocatable :: u(:,:) select case(k) case(5) ! N=31 a(1:16) = (/ -0.9999093846951440_DP, -0.9996882640283532_DP, -0.9990651964557858_DP, & -0.9975148564572244_DP, -0.9940555066314022_DP, -0.9870405605073769_DP, -0.9739668681956775_DP, & -0.9513679640727469_DP, -0.9148792632645746_DP, -0.8595690586898966_DP, -0.7806074389832004_DP, & -0.6742714922484359_DP, -0.5391467053879677_DP, -0.3772097381640342_DP, -0.1943570033249354_DP, & 0.000000000000000_DP /) a(17:31) = -a(15:1:-1) w(1:16) = (/ 0.1187484378271595D-03, 0.3628302624680807D-03, 0.9678665925484080D-03, & 0.2292994107706537D-02, 0.4897085482322761D-02, 0.9548627018272046D-02, 0.1717849940927817D-01, & 0.2875403116335220D-01, 0.4505960770348939D-01, 0.6638762853525473D-01, 0.9218368027840027D-01, & 0.1207522455632202_DP, 0.1491892696359569_DP, 0.1737092867706880_DP, 0.1904186225373517_DP, & 0.1963579530037270_DP /) w(17:31) = w(15:1:-1) case(6) ! N=63 a(1:32) = (/ -0.9999538710056279_DP, & -0.9999093846951440_DP, -0.9998288220728749_DP, -0.9996882640283532_DP, -0.9994514344352746_DP, & -0.9990651964557858_DP, -0.9984542087676977_DP, -0.9975148564572244_DP, -0.9961086654375085_DP, & -0.9940555066314022_DP, -0.9911269924416988_DP, -0.9870405605073769_DP, -0.9814548266773352_DP, & -0.9739668681956775_DP, -0.9641121642235473_DP, -0.9513679640727469_DP, -0.9351608575219846_DP, & -0.9148792632645746_DP, -0.8898914027842602_DP, -0.8595690586898966_DP, -0.8233170055064023_DP, & -0.7806074389832004_DP, -0.7310180347925614_DP, -0.6742714922484359_DP, -0.6102736575006389_DP, & -0.5391467053879677_DP, -0.4612535439395857_DP, -0.3772097381640342_DP, -0.2878799327427159_DP, & -0.1943570033249354_DP, -0.9792388528783233D-01, 0.000000000000000_DP /) a(33:63) = -a(31:1:-1) w(1:32) = (/ 0.3208813366945178D-04, & 0.5937356198247694D-04, 0.1056829831739953D-03, 0.1814131240119765D-03, 0.3010281997741137D-03, & 0.4839279419160360D-03, 0.7552454743965616D-03, 0.1146484368725149D-02, 0.1695898266447088D-02, & 0.2448515649876432D-02, 0.3455715046209973D-02, 0.4774260684945725D-02, 0.6464715068326731D-02, & 0.8589154671047334D-02, 0.1120813191793253D-01, 0.1437685651080538D-01, 0.1814061891180582D-01, & 0.2252955457639992D-01, 0.2755295281779945D-01, 0.3319344700302310D-01, 0.3940157301229943D-01, & 0.4609133016763458D-01, 0.5313749207959197D-01, 0.6037545476519407D-01, 0.6760433485730893D-01, & 0.7459380948441433D-01, 0.8109481135285784D-01, 0.8685368240400047D-01, 0.9162881336028068D-01, & 0.9520825784921207D-01, 0.9742643064229337D-01, 0.9817789022528611D-01 /) w(33:63) = w(31:1:-1) case(7) ! N=127 a(1:64) = (/ -0.9999675930679435_DP, & -0.9999538710056279_DP, -0.9999350199250824_DP, -0.9999093846951440_DP, -0.9998748650487803_DP, & -0.9998288220728749_DP, -0.9997679715995609_DP, -0.9996882640283532_DP, -0.9995847503515176_DP, & -0.9994514344352746_DP, -0.9992811119217919_DP, -0.9990651964557858_DP, -0.9987935342988059_DP, & -0.9984542087676977_DP, -0.9980333363154338_DP, -0.9975148564572244_DP, -0.9968803181281919_DP, & -0.9961086654375085_DP, -0.9951760261553273_DP, -0.9940555066314022_DP, -0.9927169971968273_DP, & -0.9911269924416988_DP, -0.9892484310901339_DP, -0.9870405605073769_DP, -0.9844588311674308_DP, & -0.9814548266773352_DP, -0.9779762351866650_DP, -0.9739668681956775_DP, -0.9693667328969173_DP, & -0.9641121642235473_DP, -0.9581360227102137_DP, -0.9513679640727469_DP, -0.9437347860527572_DP, & -0.9351608575219846_DP, -0.9255686340686127_DP, -0.9148792632645746_DP, -0.9030132815135739_DP, & -0.8898914027842602_DP, -0.8754353976304087_DP, -0.8595690586898966_DP, -0.8422192463507568_DP, & -0.8233170055064023_DP, -0.8027987413432413_DP, -0.7806074389832004_DP, -0.7566939086337300_DP, & -0.7310180347925614_DP, -0.7035500051471421_DP, -0.6742714922484359_DP, -0.6431767589852047_DP, & -0.6102736575006389_DP, -0.5755844906351516_DP, -0.5391467053879677_DP, -0.5010133893793091_DP, & -0.4612535439395857_DP, -0.4199521112784471_DP, -0.3772097381640342_DP, -0.3331422645776381_DP, & -0.2878799327427159_DP, -0.2415663195388837_DP, -0.1943570033249354_DP, -0.1464179842905879_DP, & -0.9792388528783233D-01, -0.4905596730507789D-01, 0.000000000000000_DP /) a(65:127) = -a(63:1:-1) w(1:64) = (/ 0.1161490439541570D-04, & 0.1604398885852894D-04, 0.2193370025503293D-04, 0.2968663670972276D-04, 0.3979256936961218D-04, & 0.5284123477065571D-04, 0.6953570627161861D-04, 0.9070612116060810D-04, 0.1173235010440466D-03, & 0.1505133683694674D-03, 0.1915688324651430D-03, 0.2419627949824722D-03, 0.3033589073744850D-03, & 0.3776209019039487D-03, 0.4668199121535836D-03, 0.5732393983329520D-03, 0.6993772958933317D-03, & 0.8479450120833562D-03, 0.1021862903213154D-02, 0.1224251874890123D-02, 0.1458420758058122D-02, & 0.1727849125498576D-02, 0.2036165227271721D-02, 0.2387118740717046D-02, 0.2784548053021644D-02, & 0.3232341824496342D-02, 0.3734394621555303D-02, 0.4294556463333869D-02, 0.4916576198477945D-02, & 0.5604038722495898D-02, 0.6360296164565423D-02, 0.7188393318730888D-02, 0.8090987770112945D-02, & 0.9070265373050116D-02, 0.1012785197447841D-01, 0.1126472253995556D-01, 0.1248110912697219D-01, & 0.1377640945347946D-01, 0.1514909811907401D-01, 0.1659664283933321D-01, 0.1811542833488304D-01, & 0.1970069075783202D-01, 0.2134646571905346D-01, 0.2304555307896604D-01, 0.2478950166021179D-01, & 0.2656861691232506D-01, 0.2837199428819532D-01, 0.3018758066641841D-01, 0.3200226556544212D-01, & 0.3380200314583941D-01, 0.3557196509688968D-01, 0.3729672347452612D-01, 0.3896046143248236D-01, & 0.4054720861090164D-01, 0.4204109677008743D-01, 0.4342663014205159D-01, 0.4468896398360245D-01, & 0.4581418401632107D-01, 0.4678957889054890D-01, 0.4760389756250997D-01, 0.4824758356181584D-01, & 0.4871297856875025D-01, 0.4899448851409385D-01, 0.4908870653415225D-01 /) w(65:127) = w(63:1:-1) case(8) ! N=255 a(1:128) = (/ -0.9999729464252323_DP, -0.9999675930679435_DP, & -0.9999612848078566_DP, -0.9999538710056279_DP, -0.9999451806144587_DP, -0.9999350199250824_DP, & -0.9999231701292893_DP, -0.9999093846951440_DP, -0.9998933865475925_DP, -0.9998748650487803_DP, & -0.9998534727731114_DP, -0.9998288220728749_DP, -0.9998004814311384_DP, -0.9997679715995609_DP, & -0.9997307615198084_DP, -0.9996882640283532_DP, -0.9996398313456004_DP, -0.9995847503515176_DP, & -0.9995222376512172_DP, -0.9994514344352746_DP, -0.9993714011409377_DP, -0.9992811119217919_DP, & -0.9991794489348860_DP, -0.9990651964557858_DP, -0.9989370348335121_DP, -0.9987935342988059_DP, & -0.9986331486406774_DP, -0.9984542087676977_DP, -0.9982549161719962_DP, -0.9980333363154338_DP, & -0.9977873919589065_DP, -0.9975148564572244_DP, -0.9972133470434686_DP, -0.9968803181281919_DP, & -0.9965130546402537_DP, -0.9961086654375085_DP, -0.9956640768169531_DP, -0.9951760261553273_DP, & -0.9946410557125112_DP, -0.9940555066314022_DP, -0.9934155131692640_DP, -0.9927169971968273_DP, & -0.9919556630026776_DP, -0.9911269924416988_DP, -0.9902262404675277_DP, -0.9892484310901339_DP, & -0.9881883538007427_DP, -0.9870405605073769_DP, -0.9857993630252835_DP, -0.9844588311674308_DP, & -0.9830127914811011_DP, -0.9814548266773352_DP, -0.9797782758006157_DP, -0.9779762351866650_DP, & -0.9760415602565767_DP, -0.9739668681956775_DP, -0.9717445415654873_DP, -0.9693667328969173_DP, & -0.9668253703123558_DP, -0.9641121642235473_DP, -0.9612186151511164_DP, -0.9581360227102137_DP, & -0.9548554958050227_DP, -0.9513679640727469_DP, -0.9476641906151531_DP, -0.9437347860527572_DP, & -0.9395702239332747_DP, -0.9351608575219846_DP, -0.9304969379971534_DP, -0.9255686340686127_DP, & -0.9203660530319528_DP, -0.9148792632645746_DP, -0.9090983181630204_DP, -0.9030132815135739_DP, & -0.8966142542800760_DP, -0.8898914027842602_DP, -0.8828349882446689_DP, -0.8754353976304087_DP, & -0.8676831757756460_DP, -0.8595690586898966_DP, -0.8510840079878488_DP, -0.8422192463507568_DP, & -0.8329662939194109_DP, -0.8233170055064023_DP, -0.8132636085029739_DP, -0.8027987413432413_DP, & -0.7919154923761421_DP, -0.7806074389832004_DP, -0.7688686867682466_DP, -0.7566939086337300_DP, & -0.7440783835473473_DP, -0.7310180347925614_DP, -0.7175094674873241_DP, -0.7035500051471421_DP, & -0.6891377250616677_DP, -0.6742714922484359_DP, -0.6589509917433501_DP, -0.6431767589852047_DP, & -0.6269502080510430_DP, -0.6102736575006389_DP, -0.5931503535919532_DP, -0.5755844906351516_DP, & -0.5575812282607782_DP, -0.5391467053879677_DP, -0.5202880506912302_DP, -0.5010133893793091_DP, & -0.4813318461169050_DP, -0.4612535439395857_DP, -0.4407895990339009_DP, -0.4199521112784471_DP, & -0.3987541504672377_DP, -0.3772097381640342_DP, -0.3553338251650745_DP, -0.3331422645776381_DP, & -0.3106517805528459_DP, -0.2878799327427159_DP, -0.2648450765834480_DP, -0.2415663195388837_DP, & -0.2180634734697120_DP, -0.1943570033249354_DP, -0.1704679723820105_DP, -0.1464179842905879_DP, & -0.1222291222015577_DP, -0.9792388528783233D-01, -0.7352512298567129D-01,-0.4905596730507789D-01, & -0.2453976357464916D-01, 0.000000000000000_DP /) a(129:255) = -a(127:1:-1) w(1:128) = (/ 0.4921560307232998D-05, 0.5807439097198490D-05, & 0.6834464736551962D-05, 0.8021976333168366D-05, 0.9391479729414297D-05, 0.1096682538838511D-04, & 0.1277439336593531D-04, 0.1484328487115382D-04, 0.1720551983262035D-04, 0.1989623980256602D-04, & 0.2295391544676706D-04, 0.2642055778543214D-04, 0.3034193227027582D-04, 0.3476777470613483D-04, & 0.3975200795256068D-04, 0.4535295827241165D-04, 0.5163357013109941D-04, 0.5866161819225363D-04, & 0.6650991520349031D-04, 0.7525651441994497D-04, 0.8498490517338812D-04, 0.9578420016111737D-04, & 0.1077493130013294D-03, 0.1209811245801770D-03, 0.1355866366999198D-03, 0.1516791115271980D-03, & 0.1693781953350486D-03, 0.1888100250314509D-03, 0.2101073159703949D-03, 0.2334094295482540D-03, & 0.2588624190980828D-03, 0.2866190526068620D-03, 0.3168388107952311D-03, 0.3496878591154798D-03, & 0.3853389922411368D-03, 0.4239715496402142D-03, 0.4657713008437973D-03, 0.5109302990422449D-03, & 0.5596467016628039D-03, 0.6121245566051857D-03, 0.6685735528359555D-03, 0.7292087340689185D-03, & 0.7942501742877347D-03, 0.8639226138995281D-03, 0.9384550553452169D-03, 0.1018080317034812D-02, & 0.1103034544525193D-02, 0.1193556677915233D-02, 0.1289887874500138D-02, 0.1392270885805003D-02, & 0.1500949388208539D-02, 0.1616167266473494D-02, 0.1738167849622113D-02, 0.1867193098735087D-02, & 0.2003482746412242D-02, 0.2147273387814999D-02, 0.2298797523415569D-02, 0.2458282553807834D-02, & 0.2625949727191397D-02, 0.2802013040424356D-02, 0.2986678094853407D-02, 0.3180140908472528D-02, & 0.3382586686334190D-02, 0.3594188551540296D-02, 0.3815106239573258D-02, 0.4045484759190482D-02, & 0.4285453023596742D-02, 0.4535122456126910D-02, 0.4794585575214253D-02, 0.5063914563983784D-02, & 0.5343159830393102D-02, 0.5632348564440302D-02, 0.5931483299565078D-02, 0.6240540485980043D-02, & 0.6559469084277679D-02, 0.6888189188257894D-02, 0.7226590686503984D-02, 0.7574531972792687D-02, & 0.7931838715948011D-02, 0.8298302700229281D-02, 0.8673680747771299D-02, 0.9057693734958360D-02, & 0.9450025714903295D-02, 0.9850323158407353D-02, 0.1025819432588503D-01, 0.1067320878274018D-01, & 0.1109489707056434D-01, 0.1152275054628725D-01, 0.1195622140103138D-01, 0.1239472286990171D-01, & 0.1283762964327005D-01, 0.1328427848928526D-01, 0.1373396909635332D-01, 0.1418596514318194D-01, & 0.1463949560267392D-01, 0.1509375628448443D-01, 0.1554791161943442D-01, 0.1600109668720402D-01, & 0.1645241948682547D-01, 0.1690096344746955D-01, 0.1734579017488673D-01, 0.1778594242664566D-01, & 0.1822044730702644D-01, 0.1864831967010019D-01, 0.1906856571718493D-01, 0.1948018677253986D-01, & 0.1988218321887574D-01, 0.2027355857204934D-01, 0.2065332367220688D-01, 0.2102050096667823D-01, & 0.2137412885813173D-01, 0.2171326608991022D-01, 0.2203699613911310D-01, 0.2234443158689293D-01, & 0.2263471843462563D-01, 0.2290704033411088D-01, 0.2316062269978247D-01, 0.2339473667107199D-01, & 0.2360870289358053D-01, 0.2380189508857838D-01, 0.2397374338156889D-01, 0.2412373736221462D-01, & 0.2425142884981818D-01, 0.2435643434076045D-01, 0.2443843711680134D-01, 0.2449718899591576D-01, & 0.2453251171033698D-01, 0.2454429789967597D-01 /) w(129:255) = w(127:1:-1) case default !! compute weights not saved as constants above N = 2**k-1 r = (N-1)/2 h = 4.0_DP/2**k allocate(u(2,N)) u(1,1:N) = PIOV2*cosh(h*(/(i,i=-r,r)/)) u(2,1:N) = PIOV2*sinh(h*(/(i,i=-r,r)/)) a(1:N) = tanh(u(2,1:N)) w(1:N) = u(1,1:N)/cosh(u(2,:))**2 w(1:N) = 2.0_DP*w(1:N)/sum(w(1:N)) end select ! map the -1<=x<=1 interval onto 0<=a<=s a = (a + 1.0_DP)*s/2.0_DP end subroutine tanh_sinh_setup !! ################################################### subroutine gauss_lobatto_setup(ord,abcissa,weight) implicit none 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 print *, 'unsupported Gauss-Lobatto integration order selected' stop 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 : RONE, RZERO implicit none integer, parameter :: MINTERMS = 4 real(DP), dimension(1:), intent(in) :: series real(DP) :: denom, accsum integer :: np, i, j, m #ifdef DEBUG integer :: mm #endif 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(*,'(A)',advance='no') 'not enough Wynn-Epsilon series to accelerate' accsum = 1.0/0.0 ! make it clear answer is bogus goto 777 else write(*,'(A,I3,A)',advance='no') '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,A)',advance='no') 'epsilon cancel ',m,j,' ' #ifdef DEBUG write(*,*) 'ep' do mm=1,np write(*,'(ES11.3E3,1X)') ep(mm,:) end do write(*,*) 'accsum' write(*,*) accsum #endif 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 ! modified from numerical recipes f90 (section 3.1) SUBROUTINE polint(xa,ya,x,y,dy) ! xa and ya are given x and y locations to fit an nth degree polynomial ! through. results is a value y at given location x, with error estimate dy use constants, only : DP IMPLICIT NONE REAL(DP), DIMENSION(:), INTENT(IN) :: xa,ya REAL(DP), INTENT(IN) :: x REAL(DP), INTENT(OUT) :: y,dy INTEGER :: m,n,ns REAL(DP), DIMENSION(size(xa)) :: c,d,den,ho integer, dimension(1) :: tmp n=size(xa) c=ya d=ya ho=xa-x tmp = minloc(abs(x-xa)) ns = tmp(1) y=ya(ns) ns=ns-1 do m=1,n-1 den(1:n-m)=ho(1:n-m)-ho(1+m:n) if (any(den(1:n-m) == 0.0)) then print *, 'polint: calculation failure' stop end if den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m) d(1:n-m)=ho(1+m:n)*den(1:n-m) c(1:n-m)=ho(1:n-m)*den(1:n-m) if (2*ns < n-m) then dy=c(ns+1) else dy=d(ns) ns=ns-1 end if y=y+dy end do END SUBROUTINE polint end program Driver