istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

always use generic math functions

Browse Source 
This replaces double precision-specific function in order to allow
building GRAY with other real kinds, for example single or quadruple
precision.
This commit is contained in:
Michele Guerini Rocco 2023-09-20 16:03:08 +02:00
parent dcc832ce65
commit 8bc5ac0064
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
5 changed files with 99 additions and 99 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/dispersion.f90
View File

@ -1467,7 +1467,7 @@ pure subroutine zetac (xi, yi, zr, zi, iflag)
! accuracy ! accuracy
! !
qrho = (1-0.85_wp_*y)*sqrt(qrho) qrho = (1-0.85_wp_*y)*sqrt(qrho)
n = idnint(6 + 72*qrho) n = nint(6 + 72*qrho)
j = 2*n+1 j = 2*n+1
xsum = 1.0_wp_/j xsum = 1.0_wp_/j
ysum = 0.0_wp_ ysum = 0.0_wp_
@ -1503,13 +1503,13 @@ pure subroutine zetac (xi, yi, zr, zi, iflag)
h = 0.0_wp_ h = 0.0_wp_
kapn = 0 kapn = 0
qrho = sqrt(qrho) qrho = sqrt(qrho)
nu = idint(3 + (1442/(26*qrho+77))) nu = int(3 + (1442/(26*qrho+77)))
else else
qrho = (1-y)*sqrt(1-qrho) qrho = (1-y)*sqrt(1-qrho)
h = 1.88_wp_*qrho h = 1.88_wp_*qrho
h2 = 2*h h2 = 2*h
kapn = idnint(7 + 34*qrho) kapn = nint(7 + 34*qrho)
nu = idnint(16 + 26*qrho) nu = nint(16 + 26*qrho)
endif endif
if (h>0.0_wp_) qlambda = h2**kapn if (h>0.0_wp_) qlambda = h2**kapn
rx = 0.0_wp_ rx = 0.0_wp_

16
src/eccd.f90
View File

@ -210,13 +210,13 @@ contains
apar(2) = anpl apar(2) = anpl
apar(3) = amu apar(3) = amu
apar(4) = anprre apar(4) = anprre
apar(5) = dble(e(1)) apar(5) = real(e(1))
apar(6) = dimag(e(1)) apar(6) = imag(e(1))
apar(7) = dble(e(2)) apar(7) = real(e(2))
apar(8) = dimag(e(2)) apar(8) = imag(e(2))
apar(9) = dble(e(3)) apar(9) = real(e(3))
apar(10) = dimag(e(3)) apar(10) = imag(e(3))
apar(11) = dble(ithn) apar(11) = real(ithn)
npar=size(apar) npar=size(apar)
apar(nfpp+1:npar) = eccdpar apar(nfpp+1:npar) = eccdpar
@ -251,7 +251,7 @@ contains
if(duu.le.dumin) cycle if(duu.le.dumin) cycle
apar(12) = dble(nhn) apar(12) = real(nhn, wp_)
apar(13) = ygn apar(13) = ygn
call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, & call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, &

48
src/gray_core.f90
View File

@ -735,10 +735,10 @@ contains
sss = sk - ui*sw sss = sk - ui*sw
qi1 = half*(sss + ddd) qi1 = half*(sss + ddd)
qi2 = half*(sss - ddd) qi2 = half*(sss - ddd)
rci1 = dble(qi1) rci1 = real(qi1)
rci2 = dble(qi2) rci2 = real(qi2)
ww1 = -dimag(qi1) ww1 = -imag(qi1)
ww2 = -dimag(qi2) ww2 = -imag(qi2)
else else
rci1 = rcicsi rci1 = rcicsi
rci2 = rcieta rci2 = rcieta
@ -757,12 +757,12 @@ contains
qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2
qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2
qqxy = -(qi1 - qi2)*sin(phic)*cos(phic) qqxy = -(qi1 - qi2)*sin(phic)*cos(phic)
wwxx = -dimag(qqxx) wwxx = -imag(qqxx)
wwyy = -dimag(qqyy) wwyy = -imag(qqyy)
wwxy = -dimag(qqxy) wwxy = -imag(qqxy)
rcixx = dble(qqxx) rcixx = real(qqxx)
rciyy = dble(qqyy) rciyy = real(qqyy)
rcixy = dble(qqxy) rcixy = real(qqxy)
dqi1 = -qi1**2 dqi1 = -qi1**2
dqi2 = -qi2**2 dqi2 = -qi2**2
@ -775,27 +775,27 @@ contains
d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2 d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2
d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic) d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic)
dwwxx = -dimag(dqqxx) dwwxx = -imag(dqqxx)
dwwyy = -dimag(dqqyy) dwwyy = -imag(dqqyy)
dwwxy = -dimag(dqqxy) dwwxy = -imag(dqqxy)
d2wwxx = -dimag(d2qqxx) d2wwxx = -imag(d2qqxx)
d2wwyy = -dimag(d2qqyy) d2wwyy = -imag(d2qqyy)
d2wwxy = -dimag(d2qqxy) d2wwxy = -imag(d2qqxy)
drcixx = dble(dqqxx) drcixx = real(dqqxx)
drciyy = dble(dqqyy) drciyy = real(dqqyy)
drcixy = dble(dqqxy) drcixy = real(dqqxy)
if(nrayr > 1) then if(nrayr > 1) then
dr = rwmax/dble(nrayr-1) dr = rwmax/(nrayr-1)
else else
dr = one dr = one
end if end if
ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/dble(nrayr-1) ddfu = two*dr**2/ak0 ! twodr2 = 2*dr**2 = 2*rwmax/real(nrayr-1)
do j = 1, nrayr do j = 1, nrayr
uj(j) = dble(j-1) uj(j) = real(j-1, wp_)
end do end do
da=2*pi/dble(nrayth) da=2*pi/nrayth
do k=1,nrayth do k=1,nrayth
alfak = (k-1)*da alfak = (k-1)*da
sna(k) = sin(alfak) sna(k) = sin(alfak)
@ -2280,7 +2280,7 @@ bb: do
! compute Btot=Bres/n with n=1,5 ! compute Btot=Bres/n with n=1,5
write (ubres, *) '#i Btot R z' write (ubres, *) '#i Btot R z'
do n = 1, 5 do n = 1, 5
bbb = bres/dble(n) bbb = bres/n
if (bbb >= btmn .and. bbb <= btmx) then if (bbb >= btmn .and. bbb <= btmx) then
nconts = size(ncpts) nconts = size(ncpts)
nctot = size(rrcb) nctot = size(rrcb)

6
src/polarization.f90
View File

@ -5,15 +5,15 @@ module polarization
contains contains
subroutine stokes_ce(ext,eyt,qq,uu,vv) subroutine stokes_ce(ext,eyt,qq,uu,vv)
use const_and_precisions, only : wp_,two use const_and_precisions, only : wp_
implicit none implicit none
! arguments ! arguments
complex(wp_), intent(in) :: ext,eyt complex(wp_), intent(in) :: ext,eyt
real(wp_), intent(out) :: qq,uu,vv real(wp_), intent(out) :: qq,uu,vv
qq = abs(ext)**2 - abs(eyt)**2 qq = abs(ext)**2 - abs(eyt)**2
uu = two* dble(ext*dconjg(eyt)) uu = 2*real(ext*conjg(eyt))
vv = two*dimag(ext*dconjg(eyt)) vv = 2*imag(ext*conjg(eyt))
end subroutine stokes_ce end subroutine stokes_ce

120
src/quadpack.f90
View File

@ -357,7 +357,7 @@ contains
real(wp_), external :: f real(wp_), external :: f
! !
real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, & real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, &
defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & defabs,defab1,defab2,dres,erlarg,erlast,errbnd,errmax, &
error1,error2,erro12,errsum,ertest,resabs,reseps,small error1,error2,erro12,errsum,ertest,resabs,reseps,small
real(wp_) :: res3la(3),rlist2(52) real(wp_) :: res3la(3),rlist2(52)
integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@ -430,7 +430,7 @@ contains
blist(1) = b blist(1) = b
rlist(1) = 0.0e+00_wp_ rlist(1) = 0.0e+00_wp_
elist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_
if(epsabs<=0.0e+00_wp_.and.epsrel<dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_)) then if(epsabs<=0.0e+00_wp_.and.epsrel<max(0.5e+02_wp_*epmach,0.5e-28_wp_)) then
ier = 6 ier = 6
return return
end if end if
@ -444,7 +444,7 @@ contains
! test on accuracy. ! test on accuracy.
! !
dres = abs(result) dres = abs(result)
errbnd = dmax1(epsabs,epsrel*dres) errbnd = max(epsabs,epsrel*dres)
last = 1 last = 1
rlist(1) = result rlist(1) = result
elist(1) = abserr elist(1) = abserr
@ -508,7 +508,7 @@ contains
end if end if
rlist(maxerr) = area1 rlist(maxerr) = area1
rlist(last) = area2 rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*abs(area)) errbnd = max(epsabs,epsrel*abs(area))
! !
! test for roundoff error and eventually set error flag. ! test for roundoff error and eventually set error flag.
! !
@ -523,7 +523,7 @@ contains
! set error flag in the case of bad integrand behaviour ! set error flag in the case of bad integrand behaviour
! at a point of the integration range. ! at a point of the integration range.
! !
if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & if(max(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* &
(abs(a2)+0.1e+04_wp_*uflow)) ier = 4 (abs(a2)+0.1e+04_wp_*uflow)) ier = 4
! !
! append the newly-created intervals to the list. ! append the newly-created intervals to the list.
@ -602,7 +602,7 @@ contains
abserr = abseps abserr = abseps
result = reseps result = reseps
correc = erlarg correc = erlarg
ertest = dmax1(epsabs,epsrel*abs(reseps)) ertest = max(epsabs,epsrel*abs(reseps))
! ***jump out of do-loop ! ***jump out of do-loop
if(abserr<=ertest) exit if(abserr<=ertest) exit
! !
@ -637,7 +637,7 @@ contains
! test on divergence. ! test on divergence.
! !
110 continue 110 continue
if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & if(ksgn==(-1).and.max(abs(result),abs(area))<= &
defabs*0.1e-01_wp_) go to 130 defabs*0.1e-01_wp_) go to 130
if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
.or.errsum>abs(area)) ier = 6 .or.errsum>abs(area)) ier = 6
@ -713,7 +713,7 @@ contains
real(wp_), dimension(52), intent(inout) :: epstab real(wp_), dimension(52), intent(inout) :: epstab
real(wp_), dimension(3), intent(inout) :: res3la real(wp_), dimension(3), intent(inout) :: res3la
integer, intent(inout) :: n,nres integer, intent(inout) :: n,nres
real(wp_) :: abs,delta1,delta2,delta3,dmax1,epsinf,error, & real(wp_) :: abs,delta1,delta2,delta3,epsinf,error, &
err1,err2,err3,e0,e1,e1abs,e2,e3,res,ss,tol1,tol2,tol3 err1,err2,err3,e0,e1,e1abs,e2,e3,res,ss,tol1,tol2,tol3
integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,newelm,num integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,newelm,num
! !
@ -762,10 +762,10 @@ contains
e1abs = abs(e1) e1abs = abs(e1)
delta2 = e2-e1 delta2 = e2-e1
err2 = abs(delta2) err2 = abs(delta2)
tol2 = dmax1(abs(e2),e1abs)*epmach tol2 = max(abs(e2),e1abs)*epmach
delta3 = e1-e0 delta3 = e1-e0
err3 = abs(delta3) err3 = abs(delta3)
tol3 = dmax1(e1abs,abs(e0))*epmach tol3 = max(e1abs,abs(e0))*epmach
if(err2<=tol2.and.err3<=tol3) then if(err2<=tol2.and.err3<=tol3) then
! !
! if e0, e1 and e2 are equal to within machine ! if e0, e1 and e2 are equal to within machine
@ -782,7 +782,7 @@ contains
epstab(k1) = e1 epstab(k1) = e1
delta1 = e1-e3 delta1 = e1-e3
err1 = abs(delta1) err1 = abs(delta1)
tol1 = dmax1(e1abs,abs(e3))*epmach tol1 = max(e1abs,abs(e3))*epmach
! !
! if two elements are very close to each other, omit ! if two elements are very close to each other, omit
! a part of the table by adjusting the value of n ! a part of the table by adjusting the value of n
@ -847,7 +847,7 @@ contains
res3la(3) = result res3la(3) = result
end if end if
100 continue 100 continue
abserr = dmax1(abserr,0.5e+01_wp_*epmach*abs(result)) abserr = max(abserr,0.5e+01_wp_*epmach*abs(result))
end subroutine dqelg end subroutine dqelg
subroutine dqk21(f,a,b,result,abserr,resabs,resasc) subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
@ -907,7 +907,7 @@ contains
real(wp_), intent(in) :: a,b real(wp_), intent(in) :: a,b
real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), intent(out) :: result,abserr,resabs,resasc
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: absc,centr,abs,dhlgth,dmax1,dmin1,fc,fsum, & real(wp_) :: absc,centr,abs,dhlgth,fc,fsum, &
fval1,fval2,hlgth,resg,resk,reskh fval1,fval2,hlgth,resg,resk,reskh
real(wp_), dimension(10) :: fv1,fv2 real(wp_), dimension(10) :: fv1,fv2
integer :: j,jtw,jtwm1 integer :: j,jtw,jtwm1
@ -1027,8 +1027,8 @@ contains
resasc = resasc*dhlgth resasc = resasc*dhlgth
abserr = abs((resk-resg)*hlgth) abserr = abs((resk-resg)*hlgth)
if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) & if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) &
abserr = resasc*dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) abserr = resasc*min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
((epmach*0.5e+02_wp_)*resabs,abserr) ((epmach*0.5e+02_wp_)*resabs,abserr)
end subroutine dqk21 end subroutine dqk21
@ -1519,7 +1519,7 @@ contains
integer, dimension(limit), intent(inout) :: iord integer, dimension(limit), intent(inout) :: iord
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, & real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, &
abs,defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & abs,defabs,defab1,defab2,dres,erlarg,erlast,errbnd,errmax, &
error1,error2,erro12,errsum,ertest,resabs,reseps,small error1,error2,erro12,errsum,ertest,resabs,reseps,small
real(wp_) :: res3la(3),rlist2(52) real(wp_) :: res3la(3),rlist2(52)
integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@ -1594,7 +1594,7 @@ contains
rlist(1) = 0.0e+00_wp_ rlist(1) = 0.0e+00_wp_
elist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_
iord(1) = 0 iord(1) = 0
if(epsabs<=0.0e+00_wp_.and.epsrel<dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_)) then if(epsabs<=0.0e+00_wp_.and.epsrel<max(0.5e+02_wp_*epmach,0.5e-28_wp_)) then
ier = 6 ier = 6
return return
end if end if
@ -1620,7 +1620,7 @@ contains
elist(1) = abserr elist(1) = abserr
iord(1) = 1 iord(1) = 1
dres = abs(result) dres = abs(result)
errbnd = dmax1(epsabs,epsrel*dres) errbnd = max(epsabs,epsrel*dres)
if(abserr<=1.0e+02_wp_*epmach*defabs.and.abserr>errbnd) ier = 2 if(abserr<=1.0e+02_wp_*epmach*defabs.and.abserr>errbnd) ier = 2
if(limit==1) ier = 1 if(limit==1) ier = 1
if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. &
@ -1680,7 +1680,7 @@ contains
end if end if
rlist(maxerr) = area1 rlist(maxerr) = area1
rlist(last) = area2 rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*abs(area)) errbnd = max(epsabs,epsrel*abs(area))
! !
! test for roundoff error and eventually set error flag. ! test for roundoff error and eventually set error flag.
! !
@ -1695,7 +1695,7 @@ contains
! set error flag in the case of bad integrand behaviour ! set error flag in the case of bad integrand behaviour
! at some points of the integration range. ! at some points of the integration range.
! !
if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & if(max(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* &
(abs(a2)+0.1e+04_wp_*uflow)) ier = 4 (abs(a2)+0.1e+04_wp_*uflow)) ier = 4
! !
! append the newly-created intervals to the list. ! append the newly-created intervals to the list.
@ -1771,7 +1771,7 @@ contains
abserr = abseps abserr = abseps
result = reseps result = reseps
correc = erlarg correc = erlarg
ertest = dmax1(epsabs,epsrel*abs(reseps)) ertest = max(epsabs,epsrel*abs(reseps))
if(abserr<=ertest) exit if(abserr<=ertest) exit
end if end if
! !
@ -1805,7 +1805,7 @@ contains
! test on divergence ! test on divergence
! !
110 continue 110 continue
if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & if(ksgn==(-1).and.max(abs(result),abs(area))<= &
defabs*0.1e-01_wp_) go to 130 defabs*0.1e-01_wp_) go to 130
if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
.or.errsum>abs(area)) ier = 6 .or.errsum>abs(area)) ier = 6
@ -1902,7 +1902,7 @@ contains
integer, intent(in) :: inf integer, intent(in) :: inf
real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), intent(out) :: result,abserr,resabs,resasc
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: absc,absc1,absc2,centr,abs,dinf,dmax1,dmin1,fc,fsum, & real(wp_) :: absc,absc1,absc2,centr,abs,dinf,fc,fsum, &
fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2 fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2
real(wp_), dimension(7) :: fv1,fv2 real(wp_), dimension(7) :: fv1,fv2
integer :: j integer :: j
@ -2014,8 +2014,8 @@ contains
resabs = resabs*hlgth resabs = resabs*hlgth
abserr = abs((resk-resg)*hlgth) abserr = abs((resk-resg)*hlgth)
if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* & if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* &
dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
((epmach*0.5e+02_wp_)*resabs,abserr) ((epmach*0.5e+02_wp_)*resabs,abserr)
end subroutine dqk15i end subroutine dqk15i
@ -2440,7 +2440,7 @@ contains
oflow=>comp_huge oflow=>comp_huge
implicit none implicit none
real(wp_) :: a,abseps,abserr,alist,area,area1,area12,area2,a1, & real(wp_) :: a,abseps,abserr,alist,area,area1,area12,area2,a1, &
a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,dmax1,dmin1, & a2,b,blist,b1,b2,correc,defabs,defab1,defab2, &
dres,elist,epsabs,epsrel,erlarg,erlast,errbnd, & dres,elist,epsabs,epsrel,erlarg,erlast,errbnd, &
errmax,error1,erro12,error2,errsum,ertest,points,pts, & errmax,error1,erro12,error2,errsum,ertest,points,pts, &
resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp
@ -2525,7 +2525,7 @@ contains
level(1) = 0 level(1) = 0
npts = npts2-2 npts = npts2-2
if(npts2.lt.2.or.limit.le.npts.or.(epsabs.le.0.0_wp_.and. & if(npts2.lt.2.or.limit.le.npts.or.(epsabs.le.0.0_wp_.and. &
epsrel.lt.dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_))) ier = 6 epsrel.lt.max(0.5e+02_wp_*epmach,0.5e-28_wp_))) ier = 6
if(ier.eq.6) return if(ier.eq.6) return
! !
! if any break points are provided, sort them into an ! if any break points are provided, sort them into an
@ -2533,11 +2533,11 @@ contains
! !
sign = 1.0_wp_ sign = 1.0_wp_
if(a.gt.b) sign = -1.0_wp_ if(a.gt.b) sign = -1.0_wp_
pts(1) = dmin1(a,b) pts(1) = min(a,b)
do i = 1,npts do i = 1,npts
pts(i+1) = points(i) pts(i+1) = points(i)
end do end do
pts(npts+2) = dmax1(a,b) pts(npts+2) = max(a,b)
nint = npts+1 nint = npts+1
a1 = pts(1) a1 = pts(1)
if(npts.ne.0) then if(npts.ne.0) then
@ -2552,7 +2552,7 @@ contains
end if end if
end do end do
end do end do
if(pts(1).ne.dmin1(a,b).or.pts(nintp1).ne.dmax1(a,b)) ier = 6 if(pts(1).ne.min(a,b).or.pts(nintp1).ne.max(a,b)) ier = 6
if(ier.eq.6) return if(ier.eq.6) return
end if end if
! !
@ -2586,8 +2586,8 @@ contains
! !
last = nint last = nint
neval = 21*nint neval = 21*nint
dres = dabs(result) dres = abs(result)
errbnd = dmax1(epsabs,epsrel*dres) errbnd = max(epsabs,epsrel*dres)
if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2 if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2
if(nint.ne.1) then if(nint.ne.1) then
do i = 1,npts do i = 1,npts
@ -2659,7 +2659,7 @@ contains
errsum = errsum+erro12-errmax errsum = errsum+erro12-errmax
area = area+area12-rlist(maxerr) area = area+area12-rlist(maxerr)
if(defab1.ne.error1.and.defab2.ne.error2) then if(defab1.ne.error1.and.defab2.ne.error2) then
if(dabs(rlist(maxerr)-area12).le.0.1d-04*dabs(area12) & if(abs(rlist(maxerr)-area12).le.0.1d-04*abs(area12) &
.and.erro12.ge.0.99d+00*errmax) then .and.erro12.ge.0.99d+00*errmax) then
if(extrap) iroff2 = iroff2+1 if(extrap) iroff2 = iroff2+1
if(.not.extrap) iroff1 = iroff1+1 if(.not.extrap) iroff1 = iroff1+1
@ -2670,7 +2670,7 @@ contains
level(last) = levcur level(last) = levcur
rlist(maxerr) = area1 rlist(maxerr) = area1
rlist(last) = area2 rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*dabs(area)) errbnd = max(epsabs,epsrel*abs(area))
! !
! test for roundoff error and eventually set error flag. ! test for roundoff error and eventually set error flag.
! !
@ -2685,8 +2685,8 @@ contains
! set error flag in the case of bad integrand behaviour ! set error flag in the case of bad integrand behaviour
! at a point of the integration range ! at a point of the integration range
! !
if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* & if(max(abs(a1),abs(b2)).le.(0.1d+01+0.1d+03*epmach)* &
(dabs(a2)+0.1d+04*uflow)) ier = 4 (abs(a2)+0.1d+04*uflow)) ier = 4
! !
! append the newly-created intervals to the list. ! append the newly-created intervals to the list.
! !
@ -2758,7 +2758,7 @@ contains
abserr = abseps abserr = abseps
result = reseps result = reseps
correc = erlarg correc = erlarg
ertest = dmax1(epsabs,epsrel*dabs(reseps)) ertest = max(epsabs,epsrel*abs(reseps))
! ***jump out of do-loop ! ***jump out of do-loop
if(abserr.lt.ertest) exit if(abserr.lt.ertest) exit
end if end if
@ -2786,7 +2786,7 @@ contains
if(ierro.eq.3) abserr = abserr+correc if(ierro.eq.3) abserr = abserr+correc
if(ier.eq.0) ier = 3 if(ier.eq.0) ier = 3
if(result.ne.0.0d+00.and.area.ne.0.0d+00) then if(result.ne.0.0d+00.and.area.ne.0.0d+00) then
if(abserr/dabs(result).gt.errsum/dabs(area))go to 190 if(abserr/abs(result).gt.errsum/abs(area))go to 190
else else
if(abserr.gt.errsum)go to 190 if(abserr.gt.errsum)go to 190
if(area.eq.0.0d+00) go to 210 if(area.eq.0.0d+00) go to 210
@ -2795,9 +2795,9 @@ contains
! test on divergence. ! test on divergence.
! !
end if end if
if(ksgn.ne.(-1).or.dmax1(dabs(result),dabs(area)).gt.resabs*0.1d-01) then if(ksgn.ne.(-1).or.max(abs(result),abs(area)).gt.resabs*0.1d-01) then
if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03.or. & if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03.or. &
errsum.gt.dabs(area)) ier = 6 errsum.gt.abs(area)) ier = 6
end if end if
go to 210 go to 210
! !
@ -3180,7 +3180,7 @@ contains
real(wp_), external :: f real(wp_), external :: f
! !
real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, & real(wp_) :: abseps,area,area1,area12,area2,a1,a2,b1,b2,correc,abs, &
defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & defabs,defab1,defab2,dres,erlarg,erlast,errbnd,errmax, &
error1,error2,erro12,errsum,ertest,resabs,reseps,small error1,error2,erro12,errsum,ertest,resabs,reseps,small
real(wp_) :: res3la(3),rlist2(52) real(wp_) :: res3la(3),rlist2(52)
integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@ -3253,7 +3253,7 @@ contains
blist(1) = b blist(1) = b
rlist(1) = 0.0e+00_wp_ rlist(1) = 0.0e+00_wp_
elist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_
if(epsabs<=0.0e+00_wp_.and.epsrel<dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_)) then if(epsabs<=0.0e+00_wp_.and.epsrel<max(0.5e+02_wp_*epmach,0.5e-28_wp_)) then
ier = 6 ier = 6
return return
end if end if
@ -3267,7 +3267,7 @@ contains
! test on accuracy. ! test on accuracy.
! !
dres = abs(result) dres = abs(result)
errbnd = dmax1(epsabs,epsrel*dres) errbnd = max(epsabs,epsrel*dres)
last = 1 last = 1
rlist(1) = result rlist(1) = result
elist(1) = abserr elist(1) = abserr
@ -3331,7 +3331,7 @@ contains
end if end if
rlist(maxerr) = area1 rlist(maxerr) = area1
rlist(last) = area2 rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*abs(area)) errbnd = max(epsabs,epsrel*abs(area))
! !
! test for roundoff error and eventually set error flag. ! test for roundoff error and eventually set error flag.
! !
@ -3346,7 +3346,7 @@ contains
! set error flag in the case of bad integrand behaviour ! set error flag in the case of bad integrand behaviour
! at a point of the integration range. ! at a point of the integration range.
! !
if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & if(max(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* &
(abs(a2)+0.1e+04_wp_*uflow)) ier = 4 (abs(a2)+0.1e+04_wp_*uflow)) ier = 4
! !
! append the newly-created intervals to the list. ! append the newly-created intervals to the list.
@ -3425,7 +3425,7 @@ contains
abserr = abseps abserr = abseps
result = reseps result = reseps
correc = erlarg correc = erlarg
ertest = dmax1(epsabs,epsrel*abs(reseps)) ertest = max(epsabs,epsrel*abs(reseps))
! ***jump out of do-loop ! ***jump out of do-loop
if(abserr<=ertest) exit if(abserr<=ertest) exit
! !
@ -3460,7 +3460,7 @@ contains
! test on divergence. ! test on divergence.
! !
110 continue 110 continue
if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & if(ksgn==(-1).and.max(abs(result),abs(area))<= &
defabs*0.1e-01_wp_) go to 130 defabs*0.1e-01_wp_) go to 130
if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
.or.errsum>abs(area)) ier = 6 .or.errsum>abs(area)) ier = 6
@ -3543,7 +3543,7 @@ contains
real(wp_), dimension(np), intent(in) :: apar real(wp_), dimension(np), intent(in) :: apar
real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), intent(out) :: result,abserr,resabs,resasc
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: absc,centr,abs,dhlgth,dmax1,dmin1,fc,fsum, & real(wp_) :: absc,centr,abs,dhlgth,fc,fsum, &
fval1,fval2,hlgth,resg,resk,reskh fval1,fval2,hlgth,resg,resk,reskh
real(wp_), dimension(10) :: fv1,fv2 real(wp_), dimension(10) :: fv1,fv2
integer :: j,jtw,jtwm1 integer :: j,jtw,jtwm1
@ -3663,8 +3663,8 @@ contains
resasc = resasc*dhlgth resasc = resasc*dhlgth
abserr = abs((resk-resg)*hlgth) abserr = abs((resk-resg)*hlgth)
if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) & if(resasc/=0.0e+00_wp_.and.abserr/=0.0e+00_wp_) &
abserr = resasc*dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) abserr = resasc*min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
((epmach*0.5e+02_wp_)*resabs,abserr) ((epmach*0.5e+02_wp_)*resabs,abserr)
end subroutine dqk21mv end subroutine dqk21mv
@ -4028,7 +4028,7 @@ contains
integer, dimension(limit), intent(inout) :: iord integer, dimension(limit), intent(inout) :: iord
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, & real(wp_) :: abseps,area,area1,area12,area2,a1,a2,boun,b1,b2,correc, &
abs,defabs,defab1,defab2,dmax1,dres,erlarg,erlast,errbnd,errmax, & abs,defabs,defab1,defab2,dres,erlarg,erlast,errbnd,errmax, &
error1,error2,erro12,errsum,ertest,resabs,reseps,small error1,error2,erro12,errsum,ertest,resabs,reseps,small
real(wp_) :: res3la(3),rlist2(52) real(wp_) :: res3la(3),rlist2(52)
integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, & integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@ -4103,7 +4103,7 @@ contains
rlist(1) = 0.0e+00_wp_ rlist(1) = 0.0e+00_wp_
elist(1) = 0.0e+00_wp_ elist(1) = 0.0e+00_wp_
iord(1) = 0 iord(1) = 0
if(epsabs<=0.0e+00_wp_.and.epsrel<dmax1(0.5e+02_wp_*epmach,0.5e-28_wp_)) then if(epsabs<=0.0e+00_wp_.and.epsrel<max(0.5e+02_wp_*epmach,0.5e-28_wp_)) then
ier = 6 ier = 6
return return
end if end if
@ -4129,7 +4129,7 @@ contains
elist(1) = abserr elist(1) = abserr
iord(1) = 1 iord(1) = 1
dres = abs(result) dres = abs(result)
errbnd = dmax1(epsabs,epsrel*dres) errbnd = max(epsabs,epsrel*dres)
if(abserr<=1.0e+02_wp_*epmach*defabs.and.abserr>errbnd) ier = 2 if(abserr<=1.0e+02_wp_*epmach*defabs.and.abserr>errbnd) ier = 2
if(limit==1) ier = 1 if(limit==1) ier = 1
if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. & if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. &
@ -4189,7 +4189,7 @@ contains
end if end if
rlist(maxerr) = area1 rlist(maxerr) = area1
rlist(last) = area2 rlist(last) = area2
errbnd = dmax1(epsabs,epsrel*abs(area)) errbnd = max(epsabs,epsrel*abs(area))
! !
! test for roundoff error and eventually set error flag. ! test for roundoff error and eventually set error flag.
! !
@ -4204,7 +4204,7 @@ contains
! set error flag in the case of bad integrand behaviour ! set error flag in the case of bad integrand behaviour
! at some points of the integration range. ! at some points of the integration range.
! !
if(dmax1(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* & if(max(abs(a1),abs(b2))<=(0.1e+01_wp_+0.1e+03_wp_*epmach)* &
(abs(a2)+0.1e+04_wp_*uflow)) ier = 4 (abs(a2)+0.1e+04_wp_*uflow)) ier = 4
! !
! append the newly-created intervals to the list. ! append the newly-created intervals to the list.
@ -4280,7 +4280,7 @@ contains
abserr = abseps abserr = abseps
result = reseps result = reseps
correc = erlarg correc = erlarg
ertest = dmax1(epsabs,epsrel*abs(reseps)) ertest = max(epsabs,epsrel*abs(reseps))
if(abserr<=ertest) exit if(abserr<=ertest) exit
end if end if
! !
@ -4314,7 +4314,7 @@ contains
! test on divergence ! test on divergence
! !
110 continue 110 continue
if(ksgn==(-1).and.dmax1(abs(result),abs(area))<= & if(ksgn==(-1).and.max(abs(result),abs(area))<= &
defabs*0.1e-01_wp_) go to 130 defabs*0.1e-01_wp_) go to 130
if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ & if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
.or.errsum>abs(area)) ier = 6 .or.errsum>abs(area)) ier = 6
@ -4416,7 +4416,7 @@ contains
real(wp_), dimension(np), intent(in) :: apar real(wp_), dimension(np), intent(in) :: apar
real(wp_), intent(out) :: result,abserr,resabs,resasc real(wp_), intent(out) :: result,abserr,resabs,resasc
real(wp_), external :: f real(wp_), external :: f
real(wp_) :: absc,absc1,absc2,centr,abs,dinf,dmax1,dmin1,fc,fsum, & real(wp_) :: absc,absc1,absc2,centr,abs,dinf,fc,fsum, &
fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2 fval1,fval2,hlgth,resg,resk,reskh,tabsc1,tabsc2
real(wp_), dimension(7) :: fv1,fv2 real(wp_), dimension(7) :: fv1,fv2
integer :: j integer :: j
@ -4528,8 +4528,8 @@ contains
resabs = resabs*hlgth resabs = resabs*hlgth
abserr = abs((resk-resg)*hlgth) abserr = abs((resk-resg)*hlgth)
if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* & if(resasc/=0.0e+00_wp_.and.abserr/=0._wp_) abserr = resasc* &
dmin1(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_) min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 & if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
((epmach*0.5e+02_wp_)*resabs,abserr) ((epmach*0.5e+02_wp_)*resabs,abserr)
end subroutine dqk15imv end subroutine dqk15imv