! ! 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 ! module gauss_points implicit none private public :: Gauss_Quad_Setup contains !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% subroutine Gauss_Quad_Setup(intOrd,coord,weight) use constants, only : I4B, EP integer(I4B), intent(in) :: intOrd ! Gauss integration order real(EP), intent(out), dimension(:) :: coord,weight ! from Stroud & Secrest "Gaussian Quadrature Formulas", 1966: Table 1 ! compared to quad-precision output from GAULEG() subroutine ! in Antia numerical methods book (weights and abscissa agree to >= 30 places) select case (intOrd) case(1) coord(1) = 0.0_EP weight(1) = 2.0_EP case(2) coord(1) = 0.577350269189625764509148780502_EP coord(2) = -coord(1) weight(1:2) = 1.0_EP case(3) coord(1:2) = (/ 0.774596669241483377035853079956_EP, 0.0_EP /) coord(3) = -coord(1) weight(1:2) = (/ 0.55555555555555555555555555556_EP, & & 0.88888888888888888888888888889_EP /) weight(3) = weight(1) case(4) coord(1:2) = (/ 0.861136311594052575223946488893_EP, & & 0.339981043584856264802665759103_EP /) coord(3:4) = -coord(2:1:-1) weight(1:2) = (/ 0.347854845137453857373063949222_EP, & & 0.652145154862546142626936050778_EP /) weight(3:4) = weight(2:1:-1) case(5) coord(1:3) = (/ 0.906179845938663992797626878299_EP, & & 0.538469310105683091036314420700_EP, 0.0_EP /) coord(4:5) = -coord(2:1:-1) weight(1:3) = (/ 0.236926885056189087514264040720_EP, & & 0.478628670499366468041291514836_EP, & & 0.568888888888888888888888888889_EP /) weight(4:5) = weight(2:1:-1) case(6) coord(1:3) = (/ 0.932469514203152027812301554494_EP, & & 0.661209386466264513661399595020_EP, & & 0.238619186083196908630501721681_EP /) coord(4:6) = -coord(3:1:-1) weight(1:3) = (/ 0.171324492379170345040296142173_EP, & & 0.360761573048138607569833513838_EP, & & 0.467913934572691047389870343990_EP /) weight(4:6) = weight(3:1:-1) case(7) coord(1:4) = (/ 0.949107912342758524526189684048_EP, & & 0.741531185599394439863864773281_EP, & & 0.405845151377397166906606412077_EP, 0.0_EP /) coord(5:7) = -coord(3:1:-1) weight(1:4) = (/ 0.129484966168869693270611432679_EP, & & 0.279705391489276667901467771424_EP, & & 0.381830050505118944950369775489_EP, & & 0.417959183673469387755102040816_EP /) weight(5:7) = weight(3:1:-1) case(8) coord(1:4) = (/ 0.960289856497536231683560868569_EP, & & 0.796666477413626739591553936476_EP, & & 0.525532409916328985817739049189_EP, & & 0.183434642495649804939476142360_EP /) coord(5:8) = -coord(4:1:-1) weight(1:4) = (/ 0.101228536290376259152531354310_EP, & & 0.222381034453374470544355994426_EP, & & 0.313706645877887287337962201987_EP, & & 0.362683783378361982965150449277_EP /) weight(5:8) = weight(4:1:-1) case(9) coord(1:5) = (/ 0.968160239507626089835576202904_EP, & & 0.836031107326635794299429788070_EP, & & 0.613371432700590397308702039341_EP, & & 0.324253423403808929038538014643_EP, 0.0_EP /) coord(6:9) = -coord(4:1:-1) weight(1:5) = (/ 0.0812743883615744119718921581105_EP, & & 0.180648160694857404058472031243_EP, & & 0.260610696402935462318742869419_EP, & & 0.312347077040002840068630406584_EP, & & 0.330239355001259763164525069287_EP /) weight(6:9) = weight(4:1:-1) case(10) coord(1:5) = (/ 0.973906528517171720077964012084_EP, & & 0.865063366688984510732096688423_EP, & & 0.679409568299024406234327365115_EP, & & 0.433395394129247190799265943166_EP, & & 0.148874338981631210884826001130_EP /) coord(6:10) = -coord(5:1:-1) weight(1:5) = (/ 0.0666713443086881375935688098933_EP, & & 0.149451349150580593145776339658_EP, & & 0.219086362515982043995534934228_EP, & & 0.269266719309996355091226921569_EP, & & 0.295524224714752870173892994651_EP /) weight(6:10) = weight(5:1:-1) case(12) coord(1:6) = (/ 0.981560634246719250690549090149_EP, & & 0.904117256370474856678465866119_EP, & & 0.769902674194304687036893833213_EP, & & 0.587317954286617447296702418941_EP, & & 0.367831498998180193752691536644_EP, & & 0.125233408511468915472441369464_EP /) coord(7:12) = -coord(6:1:-1) weight(1:6) = (/ 0.0471753363865118271946159614850_EP, & & 0.106939325995318430960254718194_EP, & & 0.160078328543346226334652529543_EP, & & 0.203167426723065921749064455810_EP, & & 0.233492536538354808760849898925_EP, & & 0.249147045813402785000562436043_EP /) weight(7:12) = weight(6:1:-1) case(20) coord(1:10) = (/ 0.993128599185094924786122388471_EP, & & 0.963971927277913791267666131197_EP, & & 0.912234428251325905867752441203_EP, & & 0.839116971822218823394529061702_EP, & & 0.746331906460150792614305070356_EP, & & 0.636053680726515025452836696226_EP, & & 0.510867001950827098004364050955_EP, & & 0.373706088715419560672548177025_EP, & & 0.227785851141645078080496195369_EP, & & 0.0765265211334973337546404093988_EP /) coord(11:20) = -coord(10:1:-1) weight(1:10) = (/ 0.0176140071391521183118619623519_EP, & & 0.0406014298003869413310399522749_EP, & & 0.0626720483341090635695065351870_EP, & & 0.0832767415767047487247581432220_EP, & & 0.101930119817240435036750135480_EP, & & 0.118194531961518417312377377711_EP, & & 0.131688638449176626898494499748_EP, & & 0.142096109318382051329298325067_EP, & & 0.149172986472603746787828737002_EP, & & 0.152753387130725850698084331955_EP /) weight(11:20) = weight(10:1:-1) case default print *, 'unsupported Gaussian Quadrature integration order:', intOrd stop end select end subroutine Gauss_Quad_Setup end module gauss_points