! ! Copyright (c) 2008 Kristopher L. Kuhlman (kuhlman at hwr dot arizona dot edu) ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the Software is ! furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN ! THE SOFTWARE. ! ! Malama, B., K.L. Kuhlman, and W. Barrash, 2007. Semi-analytical solution ! for flow in leaky unconfined aquifer-aquitard systems, Journal of ! Hydrology, 346(1-2), 59–68. ! http://dx.doi.org/10.1016/j.jhydrol.2007.08.018 ! ! Malama, B., K.L. Kuhlman, and W. Barrash, 2008. Semi-analytical solution ! for flow in a leaky unconfined aquifer toward a partially penetrating ! pumping well, Journal of Hydrology, 356(1-2), 234–244. ! http://dx.doi.org/10.1016/j.jhydrol.2008.03.029 ! ! $Id: utility.f90,v 1.5 2007/02/19 02:21:26 kris Exp kris $ !! this file has _four_ 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) ! 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 :: 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 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,Kr,Kz,Ss real(DP), save :: Sy,zd,rd, lambda, ld,dd ! aquifer parameters (secondary - computed from above) real(DP), save :: kappa,alphaDy,bd real(DP), save, dimension(MAXNLAYERS) :: alphar,alphaz,alphaDr,alphaDz,KDz real(DP), save :: ztop, zbot, bd2,bd3 logical, save :: piezometer logical, dimension(MAXNLAYERS) :: layer ! inverse Laplace transform parameters real(DP) :: lapalpha, laptol integer :: lapM ! Gauss Quadrature local abscissa, weights, and abscissa in global coords real(EP), allocatable :: obsGQx(:), obsGQw(:), obsGQy(:) complex(EP), allocatable :: obsGQz(:,:) end module shared_data !! ###################################################################### !! ###################################################################### module complex_hyp_trig use constants, only : EP implicit none private public :: cecosh, cesinh, cecoth, cesinhcosh contains !!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !! extended-precision complex cosh pure elemental function cecosh(z) result(c) use constants, only : EHALF complex(EP), intent(in) :: z complex(EP) :: c c = EHALF*(exp(z) + exp(-z)) end function cecosh !! extended-precision complex sinh pure elemental function cesinh(z) result(s) use constants, only : EHALF complex(EP), intent(in) :: z complex(EP) :: s s = EHALF*(exp(z) - exp(-z)) end function cesinh !! extended-precision complex sinh and cosh !! for same argument pure elemental subroutine cesinhcosh(z,s,c) use constants, only : EHALF complex(EP), intent(in) :: z complex(EP), intent(out) :: s,c complex(EP) :: expplus, expmin expplus = exp(z); expmin = exp(-z) c = EHALF*(expplus + expmin) s = EHALF*(expplus - expmin) end subroutine cesinhcosh !! extended-precision complex coth pure elemental function cecoth(z) result(c) use constants, only : EONE,ETWO complex(EP), intent(in) :: z complex(EP) :: c, temp temp = exp(-ETWO*z) c = (EONE + temp)/(EONE - temp) end function cecoth end module complex_hyp_trig module invlap_fcns implicit none private public :: invlap_setup, check_nan contains !! consolidate calculation of laplace parameter vector subroutine invlap_setup(p,tee) use constants, only : DP,EP,PIEP,ETWO use shared_data, only : lapM,lapAlpha,lapTol,tsval complex(EP), intent(out), dimension(2*lapM+1) :: p real(DP), intent(out) :: tee integer :: np,i real(EP), dimension(0:2*lapM) :: run run = real((/ (i,i=0,2*lapM) /),EP) np = 2*lapM+1 tee = 2.0_DP*tsval p(1:np) = cmplx(real(lapalpha,EP) - log(real(laptol,EP))/(ETWO*tee), run*PIEP/tee, EP) ! check for potential problems if(minval(abs(p)) < lapalpha) then stop 'decrease alpha parameter for de Hoog routine' end if end subroutine invlap_setup subroutine check_nan(fp,np,M) use constants, only : EP complex(EP), intent(in), dimension(:) :: fp integer, intent(out) :: np, M logical, dimension(size(fp)) :: mask integer :: nfp, j nfp = size(fp) mask = .false. np = nfp M = (np-1)/2 where(fp /= fp) mask = .true. end where loop: do j=1,nfp if(mask(j)) then np = j exit loop end if end do loop if(np /= nfp) then if(np <= 4) then print *, 'insufficient terms for useful Laplace inversion',np else print *, 'reduced np:',nfp,'->',np end if if(mod(np,2) == 1) then !! if new np is odd M = (np-1)/2 else !! if new np is even, lose additional term M = (np-2)/2 np = np-1 end if end if end subroutine check_nan end module invlap_fcns