! ! 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 ! ! $Id: utility.f90,v 1.3 2009/03/15 03:38:31 kris Exp kris $ !! this file has _three_ modules in it !! ###################################################################### !! ###################################################################### module constants ! integers with range of 9 and 18 orders of magnitude integer, parameter :: I4B = selected_int_kind(r=9) ! std 32-bit integer ! real with range 300 orders of mag, 15 sig figs (8 on both g95 & ifort) integer, parameter :: DP = selected_real_kind (p=15,r=300) !! extended range internal variables (10 on g95, 16 on ifort) integer, parameter :: EP = selected_real_kind(r=3000) !! full quad precision integer, parameter :: QP = selected_real_kind(p=33,r=3000) ! useful? constants related to pi and ln ! calculated to precision=33 using ifort or mathematica real(DP), parameter :: PI = 3.141592653589793238462643383279503_DP real(DP), parameter :: TWOPI = 6.28318530717958647692528676655901_DP real(DP), parameter :: PIOV2 = 1.57079632679489661923132169163975_DP real(DP), parameter :: PIOV3 = 1.04719755119659774615421446109317_DP real(DP), parameter :: PIOV4 = 0.785398163397448309615660845819876_DP real(DP), parameter :: ONEOVPI = 0.318309886183790671537767526745029_DP real(DP), parameter :: LN10 = 2.30258509299404568401799145468436_DP real(DP), parameter :: LN2 = 0.693147180559945309417232121458177_DP real(DP), parameter :: LNPIOV2 = 0.451582705289454864726195229894882_DP real(DP), parameter :: EULER = 0.5772156649015328606065120900824025_DP real(DP), parameter :: PISQ = 9.86960440108935861883449099987615_DP real(DP), parameter :: SQRTPI = 1.77245385090551602729816748334115_DP real(DP), parameter :: SQRTPIOV2 = 1.25331413731550025120788264240552_DP real(DP), parameter :: INVLOG10 = 0.434294481903251827651128918916605_DP real(EP), parameter :: LN2EP = 0.693147180559945309417232121458177_EP real(QP), parameter :: LN2QP = 0.693147180559945309417232121458177_QP real(EP), parameter :: PIEP = 3.141592653589793238462643383279503_EP ! a very small number (if < SMALL, it is effectively zero, compared to one) real(DP), parameter :: SMALL = epsilon(1.0_DP) ! approx 2.22E-16 ! the largest number which can be represented accurately real(DP), parameter :: LARGE = huge(1.0_DP) ! approx 1.79E+308 ! commonly(?) used constants real(EP), parameter :: EONE = 1.0_EP real(EP), parameter :: EZERO = 0.0_EP real(EP), parameter :: EHALF = 0.5_EP real(EP), parameter :: ETWO = 2.0_EP complex(DP), parameter :: CZERO = (0.0_DP, 0.0_DP) complex(DP), parameter :: CONE = (1.0_DP, 0.0_DP) complex(DP), parameter :: CTWO = (2.0_DP, 0.0_DP) complex(DP), parameter :: EYE = (0.0_DP, 1.0_DP) real(DP), parameter :: RZERO = 0.0_DP real(DP), parameter :: RONE = 1.0_DP real(DP), parameter :: RTWO = 2.0_DP real(DP), parameter :: HALF = 0.5_DP end module constants !! ###################################################################### !! ###################################################################### module shared_data use constants, only : DP,EP,QP implicit none integer, parameter :: MAXNLAYERS = 3 ! variables to pass using modules ! parameters to change between calculations real(DP), save :: tsval ! time - passed through hankel integral to Lap xform ! aquifer parameters (primary - read from input) real(DP), save, dimension(MAXNLAYERS) :: b,sigma,gamell real(DP), save :: Kr,Kz,Ss,Sy,zd,rd ! aquifer parameters (secondary - computed from above) real(DP), save :: alphaR, alphaY, alphaD, kappa real(DP), save, dimension(1:3) :: sigmaD real(DP), save :: bd1, bd3 logical, save, dimension(MAXNLAYERS) :: layer integer, save :: stehfestN real(QP), save, allocatable :: stehfestV(:) end module shared_data !! ###################################################################### !! ###################################################################### module invlap_fcns implicit none private public :: stehfest_setup, stehfest, stehfest_coeff contains function stehfest_coeff(N) result(V) use constants, only : QP integer, intent(in) :: N real(QP), dimension(N) :: V integer :: num real :: x if (mod(N,2) == 1) then print *, 'order of stehfest coefficient must be even,',N stop end if select case(N) case(14) V = (/ 0.002777777777777777777777777777777777777778_QP,& & -6.402777777777777777777777777777777777778_QP, & & 924.05_QP, & & -34597.92777777777777777777777777777777778_QP, & & 540321.1111111111111111111111111111111111_QP, & & -4398346.366666666666666666666666666666667_QP, & & 21087591.77777777777777777777777777777778_QP, & & -63944913.04444444444444444444444444444444_QP, & & 127597579.55_QP, & & -170137188.0833333333333333333333333333333_QP,& & 150327467.0333333333333333333333333333333_QP, & & -84592161.5_QP, & & 27478884.76666666666666666666666666666667_QP, & & -3925554.966666666666666666666666666666667_QP /) case(16) V = (/ -0.0003968253968253968253968253968253968253968_QP, & & 2.133730158730158730158730158730158730159_QP, & & -551.0166666666666666666666666666666666667_QP, & & 33500.16111111111111111111111111111111111_QP, & & -812665.1111111111111111111111111111111111_QP, & & 10076183.76666666666666666666666666666667_QP, & & -73241382.97777777777777777777777777777778_QP, & & 339059632.0730158730158730158730158730159_QP, & & -1052539536.278571428571428571428571428571_QP, & & 2259013328.583333333333333333333333333333_QP, & & -3399701984.433333333333333333333333333333_QP, & & 3582450461.7_QP, & & -2591494081.366666666666666666666666666667_QP, & & 1227049828.766666666666666666666666666667_QP, & & -342734555.4285714285714285714285714285714_QP, & & 42841819.42857142857142857142857142857143_QP /) case(18) V = (/ 0.0000496031746031746031746031746031746031746_QP, & & -0.6095734126984126984126984126984126984127_QP, & & 274.5940476190476190476190476190476190476_QP, & & -26306.95674603174603174603174603174603175_QP, & & 957257.2013888888888888888888888888888889_QP, & & -17358694.84583333333333333333333333333333_QP, & & 182421222.6472222222222222222222222222222_QP, & & -1218533288.309126984126984126984126984127_QP, & & 5491680025.283035714285714285714285714286_QP, & & -17362131115.2068452380952380952380952381_QP, & & 39455096903.52738095238095238095238095238_QP, & & -65266516985.175_QP, & & 78730068328.22083333333333333333333333333_QP, & & -68556444196.12083333333333333333333333333_QP, & & 41984343475.05357142857142857142857142857_QP, & & -17160934711.83928571428571428571428571429_QP, & & 4204550039.102678571428571428571428571429_QP, & & -467172226.5669642857142857142857142857143_QP /) case default open(unit=27,file='stehfest_coefficients.dat',action='read',status='old') readV: do read(27,*) num if (num == N) then read(27,*) V exit readV elseif (num > N) then stop 'stehfestN used not supported' else read(27,*) x ! throw away this data, onto end if end do readV close(27) end select end function stehfest_coeff subroutine stehfest_setup(t,p) use constants, only : DP,EP,LN2EP use shared_data, only : stehfestN real(EP), dimension(stehfestN), intent(out) :: p real(DP), intent(in) :: t integer :: i p = LN2EP*real((/(i, i=1,stehfestN)/),EP)/t end subroutine stehfest_setup function stehfest(fp,t) result(f) use constants, only : DP,EP,LN2EP use shared_data, only : stehfestV real(EP), intent(in), dimension(:) :: fp real(DP), intent(in) :: t real(DP) :: f f = real(LN2EP/t*sum(stehfestV*fp(:)),DP) end function stehfest end module invlap_fcns