diff --git a/src/dispersion.f90 b/src/dispersion.f90
index 81ee716..a4228ad 100644
--- a/src/dispersion.f90
+++ b/src/dispersion.f90
@@ -1467,7 +1467,7 @@ pure subroutine zetac (xi, yi, zr, zi, iflag)
 !  accuracy
 !
       qrho  = (1-0.85_wp_*y)*sqrt(qrho)
-      n     = idnint(6 + 72*qrho)
+      n     = nint(6 + 72*qrho)
       j     = 2*n+1
       xsum  = 1.0_wp_/j
       ysum  = 0.0_wp_
@@ -1503,13 +1503,13 @@ pure subroutine zetac (xi, yi, zr, zi, iflag)
         h    = 0.0_wp_
         kapn = 0
         qrho = sqrt(qrho)
-        nu   = idint(3 + (1442/(26*qrho+77)))
+        nu   = int(3 + (1442/(26*qrho+77)))
       else
         qrho = (1-y)*sqrt(1-qrho)
         h    = 1.88_wp_*qrho
         h2   = 2*h
-        kapn = idnint(7  + 34*qrho)
-        nu   = idnint(16 + 26*qrho)
+        kapn = nint(7  + 34*qrho)
+        nu   = nint(16 + 26*qrho)
       endif
       if (h>0.0_wp_) qlambda = h2**kapn
       rx = 0.0_wp_
diff --git a/src/eccd.f90 b/src/eccd.f90
index b578669..d59caa2 100644
--- a/src/eccd.f90
+++ b/src/eccd.f90
@@ -210,13 +210,13 @@ contains
     apar(2) = anpl
     apar(3) = amu
     apar(4) = anprre
-    apar(5) = dble(e(1))
-    apar(6) = dimag(e(1))
-    apar(7) = dble(e(2))
-    apar(8) = dimag(e(2))
-    apar(9) = dble(e(3))
-    apar(10) = dimag(e(3))
-    apar(11) = dble(ithn)
+    apar(5) = real(e(1))
+    apar(6) = imag(e(1))
+    apar(7) = real(e(2))
+    apar(8) = imag(e(2))
+    apar(9) = real(e(3))
+    apar(10) = imag(e(3))
+    apar(11) = real(ithn)
 
     npar=size(apar)
     apar(nfpp+1:npar) = eccdpar
@@ -251,7 +251,7 @@ contains
 
       if(duu.le.dumin) cycle
 
-      apar(12) = dble(nhn)
+      apar(12) = real(nhn, wp_)
       apar(13) = ygn
 
       call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp, &
diff --git a/src/gray_core.f90 b/src/gray_core.f90
index 42365b1..b7ef3b8 100644
--- a/src/gray_core.f90
+++ b/src/gray_core.f90
@@ -735,10 +735,10 @@ contains
       sss  =  sk - ui*sw
       qi1  =  half*(sss + ddd)
       qi2  =  half*(sss - ddd)
-      rci1 =  dble(qi1)
-      rci2 =  dble(qi2)
-      ww1  = -dimag(qi1)
-      ww2  = -dimag(qi2)
+      rci1 =  real(qi1)
+      rci2 =  real(qi2)
+      ww1  = -imag(qi1)
+      ww2  = -imag(qi2)
     else
       rci1 = rcicsi
       rci2 = rcieta
@@ -757,12 +757,12 @@ contains
     qqxx  =  qi1*cos(phic)**2 + qi2*sin(phic)**2
     qqyy  =  qi1*sin(phic)**2 + qi2*cos(phic)**2
     qqxy  = -(qi1 - qi2)*sin(phic)*cos(phic)
-    wwxx  = -dimag(qqxx)
-    wwyy  = -dimag(qqyy)
-    wwxy  = -dimag(qqxy)
-    rcixx =  dble(qqxx)
-    rciyy =  dble(qqyy)
-    rcixy =  dble(qqxy)
+    wwxx  = -imag(qqxx)
+    wwyy  = -imag(qqyy)
+    wwxy  = -imag(qqxy)
+    rcixx =  real(qqxx)
+    rciyy =  real(qqyy)
+    rcixy =  real(qqxy)
 
     dqi1  =  -qi1**2
     dqi2  =  -qi2**2
@@ -775,27 +775,27 @@ contains
     d2qqyy =  d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2
     d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic)
 
-    dwwxx  = -dimag(dqqxx)
-    dwwyy  = -dimag(dqqyy)
-    dwwxy  = -dimag(dqqxy)
-    d2wwxx = -dimag(d2qqxx)
-    d2wwyy = -dimag(d2qqyy)
-    d2wwxy = -dimag(d2qqxy)
-    drcixx =  dble(dqqxx)
-    drciyy =  dble(dqqyy)
-    drcixy =  dble(dqqxy)
+    dwwxx  = -imag(dqqxx)
+    dwwyy  = -imag(dqqyy)
+    dwwxy  = -imag(dqqxy)
+    d2wwxx = -imag(d2qqxx)
+    d2wwyy = -imag(d2qqyy)
+    d2wwxy = -imag(d2qqxy)
+    drcixx =  real(dqqxx)
+    drciyy =  real(dqqyy)
+    drcixy =  real(dqqxy)
 
     if(nrayr > 1) then
-      dr = rwmax/dble(nrayr-1)
+      dr = rwmax/(nrayr-1)
     else
       dr = one
     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
-      uj(j)     = dble(j-1)
+      uj(j) = real(j-1, wp_)
     end do
 
-    da=2*pi/dble(nrayth)
+    da=2*pi/nrayth
     do k=1,nrayth
       alfak  = (k-1)*da
       sna(k) = sin(alfak)
@@ -2280,7 +2280,7 @@ bb: do
     ! compute Btot=Bres/n with n=1,5
     write (ubres, *) '#i Btot  R  z'
     do n = 1, 5
-      bbb = bres/dble(n)
+      bbb = bres/n
       if (bbb >= btmn .and. bbb <= btmx) then
         nconts = size(ncpts)
         nctot = size(rrcb)
diff --git a/src/polarization.f90 b/src/polarization.f90
index 5c28598..2ed4055 100644
--- a/src/polarization.f90
+++ b/src/polarization.f90
@@ -5,15 +5,15 @@ module polarization
 
 contains
   subroutine stokes_ce(ext,eyt,qq,uu,vv)
-    use const_and_precisions, only : wp_,two
+    use const_and_precisions, only : wp_
     implicit none
 ! arguments
     complex(wp_), intent(in) :: ext,eyt
     real(wp_), intent(out) :: qq,uu,vv
 
     qq = abs(ext)**2 - abs(eyt)**2
-    uu = two* dble(ext*dconjg(eyt))
-    vv = two*dimag(ext*dconjg(eyt))
+    uu = 2*real(ext*conjg(eyt))
+    vv = 2*imag(ext*conjg(eyt))
   end subroutine stokes_ce
 
 
diff --git a/src/quadpack.f90 b/src/quadpack.f90
index 6c3aab8..42fa7ba 100644
--- a/src/quadpack.f90
+++ b/src/quadpack.f90
@@ -357,7 +357,7 @@ contains
     real(wp_), external :: f
 !
     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
     real(wp_) :: res3la(3),rlist2(52)
     integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@@ -430,7 +430,7 @@ contains
     blist(1) = b
     rlist(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
       return
     end if
@@ -444,7 +444,7 @@ contains
 !       test on accuracy.
 !
     dres = abs(result)
-    errbnd = dmax1(epsabs,epsrel*dres)
+    errbnd = max(epsabs,epsrel*dres)
     last = 1
     rlist(1) = result
     elist(1) = abserr
@@ -508,7 +508,7 @@ contains
       end if
       rlist(maxerr) = area1
       rlist(last) = area2
-      errbnd = dmax1(epsabs,epsrel*abs(area))
+      errbnd = max(epsabs,epsrel*abs(area))
 !
 !       test for roundoff error and eventually set error flag.
 !
@@ -523,7 +523,7 @@ contains
 !       set error flag in the case of bad integrand behaviour
 !       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
 !
 !       append the newly-created intervals to the list.
@@ -602,7 +602,7 @@ contains
         abserr = abseps
         result = reseps
         correc = erlarg
-        ertest = dmax1(epsabs,epsrel*abs(reseps))
+        ertest = max(epsabs,epsrel*abs(reseps))
 ! ***jump out of do-loop
         if(abserr<=ertest) exit
 !
@@ -637,7 +637,7 @@ contains
 !       test on divergence.
 !
     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
     if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
      .or.errsum>abs(area)) ier = 6
@@ -713,7 +713,7 @@ contains
     real(wp_), dimension(52), intent(inout) :: epstab
     real(wp_), dimension(3), intent(inout) :: res3la
     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
     integer i,ib,ib2,ie,indx,k1,k2,k3,limexp,newelm,num
 !
@@ -762,10 +762,10 @@ contains
       e1abs = abs(e1)
       delta2 = e2-e1
       err2 = abs(delta2)
-      tol2 = dmax1(abs(e2),e1abs)*epmach
+      tol2 = max(abs(e2),e1abs)*epmach
       delta3 = e1-e0
       err3 = abs(delta3)
-      tol3 = dmax1(e1abs,abs(e0))*epmach
+      tol3 = max(e1abs,abs(e0))*epmach
       if(err2<=tol2.and.err3<=tol3) then
 !
 !       if e0, e1 and e2 are equal to within machine
@@ -782,7 +782,7 @@ contains
       epstab(k1) = e1
       delta1 = e1-e3
       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
 !       a part of the table by adjusting the value of n
@@ -847,7 +847,7 @@ contains
       res3la(3) = result
     end if
     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
 
   subroutine dqk21(f,a,b,result,abserr,resabs,resasc)
@@ -907,7 +907,7 @@ contains
     real(wp_), intent(in) :: a,b
     real(wp_), intent(out) :: result,abserr,resabs,resasc
     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
     real(wp_), dimension(10) :: fv1,fv2
     integer :: j,jtw,jtwm1
@@ -1027,8 +1027,8 @@ contains
     resasc = resasc*dhlgth
     abserr = abs((resk-resg)*hlgth)
     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_)
-    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 &
+      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 = max &
       ((epmach*0.5e+02_wp_)*resabs,abserr)
   end subroutine dqk21
 
@@ -1519,7 +1519,7 @@ contains
     integer, dimension(limit), intent(inout) :: iord
     real(wp_), external :: f
     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
     real(wp_) :: res3la(3),rlist2(52)
     integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@@ -1594,7 +1594,7 @@ contains
     rlist(1) = 0.0e+00_wp_
     elist(1) = 0.0e+00_wp_
     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
       return
     end if
@@ -1620,7 +1620,7 @@ contains
     elist(1) = abserr
     iord(1) = 1
     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(limit==1) ier = 1
     if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. &
@@ -1680,7 +1680,7 @@ contains
       end if
       rlist(maxerr) = area1
       rlist(last) = area2
-      errbnd = dmax1(epsabs,epsrel*abs(area))
+      errbnd = max(epsabs,epsrel*abs(area))
 !
 !       test for roundoff error and eventually set error flag.
 !
@@ -1695,7 +1695,7 @@ contains
 !       set error flag in the case of bad integrand behaviour
 !       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
 !
 !       append the newly-created intervals to the list.
@@ -1771,7 +1771,7 @@ contains
         abserr = abseps
         result = reseps
         correc = erlarg
-        ertest = dmax1(epsabs,epsrel*abs(reseps))
+        ertest = max(epsabs,epsrel*abs(reseps))
         if(abserr<=ertest) exit
       end if
 !
@@ -1805,7 +1805,7 @@ contains
 !       test on divergence
 !
     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
     if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
      .or.errsum>abs(area)) ier = 6
@@ -1902,7 +1902,7 @@ contains
     integer, intent(in) :: inf
     real(wp_), intent(out) :: result,abserr,resabs,resasc
     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
     real(wp_), dimension(7) :: fv1,fv2
     integer :: j
@@ -2014,8 +2014,8 @@ contains
     resabs = resabs*hlgth
     abserr = abs((resk-resg)*hlgth)
     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_)
-    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 &
+     min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
+    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
      ((epmach*0.5e+02_wp_)*resabs,abserr)
   end subroutine dqk15i
 
@@ -2440,7 +2440,7 @@ contains
       oflow=>comp_huge
     implicit none
     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,     &
         errmax,error1,erro12,error2,errsum,ertest,points,pts,    &
         resa,resabs,reseps,result,res3la,rlist,rlist2,sign,temp
@@ -2525,7 +2525,7 @@ contains
     level(1) = 0
     npts = npts2-2
     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 any break points are provided, sort them into an
@@ -2533,11 +2533,11 @@ contains
 !
     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
       pts(i+1) = points(i)
     end do
-    pts(npts+2) = dmax1(a,b)
+    pts(npts+2) = max(a,b)
     nint = npts+1
     a1 = pts(1)
     if(npts.ne.0) then
@@ -2552,7 +2552,7 @@ contains
           end if
         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
     end if
 !
@@ -2586,8 +2586,8 @@ contains
 !
     last = nint
     neval = 21*nint
-    dres = dabs(result)
-    errbnd = dmax1(epsabs,epsrel*dres)
+    dres = abs(result)
+    errbnd = max(epsabs,epsrel*dres)
     if(abserr.le.0.1d+03*epmach*resabs.and.abserr.gt.errbnd) ier = 2
     if(nint.ne.1) then
       do i = 1,npts
@@ -2659,7 +2659,7 @@ contains
         errsum = errsum+erro12-errmax
         area = area+area12-rlist(maxerr)
         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
             if(extrap) iroff2 = iroff2+1
             if(.not.extrap) iroff1 = iroff1+1
@@ -2670,7 +2670,7 @@ contains
         level(last) = levcur
         rlist(maxerr) = area1
         rlist(last) = area2
-        errbnd = dmax1(epsabs,epsrel*dabs(area))
+        errbnd = max(epsabs,epsrel*abs(area))
 !
 !           test for roundoff error and eventually set error flag.
 !
@@ -2685,8 +2685,8 @@ contains
 !           set error flag in the case of bad integrand behaviour
 !           at a point of the integration range
 !
-        if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)*    &
-           (dabs(a2)+0.1d+04*uflow)) ier = 4
+        if(max(abs(a1),abs(b2)).le.(0.1d+01+0.1d+03*epmach)*    &
+           (abs(a2)+0.1d+04*uflow)) ier = 4
 !
 !           append the newly-created intervals to the list.
 !
@@ -2758,7 +2758,7 @@ contains
             abserr = abseps
             result = reseps
             correc = erlarg
-            ertest = dmax1(epsabs,epsrel*dabs(reseps))
+            ertest = max(epsabs,epsrel*abs(reseps))
 ! ***jump out of do-loop
             if(abserr.lt.ertest) exit
           end if
@@ -2786,7 +2786,7 @@ contains
         if(ierro.eq.3) abserr = abserr+correc
         if(ier.eq.0) ier = 3
         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
           if(abserr.gt.errsum)go to 190
           if(area.eq.0.0d+00) go to 210
@@ -2795,9 +2795,9 @@ contains
 !           test on divergence.
 !
       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.  &
-          errsum.gt.dabs(area)) ier = 6
+          errsum.gt.abs(area)) ier = 6
       end if
       go to 210
 !
@@ -3180,7 +3180,7 @@ contains
     real(wp_), external :: f
 !
     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
     real(wp_) :: res3la(3),rlist2(52)
     integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@@ -3253,7 +3253,7 @@ contains
     blist(1) = b
     rlist(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
       return
     end if
@@ -3267,7 +3267,7 @@ contains
 !       test on accuracy.
 !
     dres = abs(result)
-    errbnd = dmax1(epsabs,epsrel*dres)
+    errbnd = max(epsabs,epsrel*dres)
     last = 1
     rlist(1) = result
     elist(1) = abserr
@@ -3331,7 +3331,7 @@ contains
       end if
       rlist(maxerr) = area1
       rlist(last) = area2
-      errbnd = dmax1(epsabs,epsrel*abs(area))
+      errbnd = max(epsabs,epsrel*abs(area))
 !
 !       test for roundoff error and eventually set error flag.
 !
@@ -3346,7 +3346,7 @@ contains
 !       set error flag in the case of bad integrand behaviour
 !       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
 !
 !       append the newly-created intervals to the list.
@@ -3425,7 +3425,7 @@ contains
         abserr = abseps
         result = reseps
         correc = erlarg
-        ertest = dmax1(epsabs,epsrel*abs(reseps))
+        ertest = max(epsabs,epsrel*abs(reseps))
 ! ***jump out of do-loop
         if(abserr<=ertest) exit
 !
@@ -3460,7 +3460,7 @@ contains
 !       test on divergence.
 !
     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
     if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
      .or.errsum>abs(area)) ier = 6
@@ -3543,7 +3543,7 @@ contains
     real(wp_), dimension(np), intent(in) :: apar
     real(wp_), intent(out) :: result,abserr,resabs,resasc
     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
     real(wp_), dimension(10) :: fv1,fv2
     integer :: j,jtw,jtwm1
@@ -3663,8 +3663,8 @@ contains
     resasc = resasc*dhlgth
     abserr = abs((resk-resg)*hlgth)
     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_)
-    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 &
+      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 = max &
       ((epmach*0.5e+02_wp_)*resabs,abserr)
   end subroutine dqk21mv
 
@@ -4028,7 +4028,7 @@ contains
     integer, dimension(limit), intent(inout) :: iord
     real(wp_), external :: f
     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
     real(wp_) :: res3la(3),rlist2(52)
     integer :: id,ierro,iroff1,iroff2,iroff3,jupbnd,k,ksgn, &
@@ -4103,7 +4103,7 @@ contains
     rlist(1) = 0.0e+00_wp_
     elist(1) = 0.0e+00_wp_
     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
       return
     end if
@@ -4129,7 +4129,7 @@ contains
     elist(1) = abserr
     iord(1) = 1
     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(limit==1) ier = 1
     if(ier/=0.or.(abserr<=errbnd.and.abserr/=resabs).or. &
@@ -4189,7 +4189,7 @@ contains
       end if
       rlist(maxerr) = area1
       rlist(last) = area2
-      errbnd = dmax1(epsabs,epsrel*abs(area))
+      errbnd = max(epsabs,epsrel*abs(area))
 !
 !       test for roundoff error and eventually set error flag.
 !
@@ -4204,7 +4204,7 @@ contains
 !       set error flag in the case of bad integrand behaviour
 !       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
 !
 !       append the newly-created intervals to the list.
@@ -4280,7 +4280,7 @@ contains
         abserr = abseps
         result = reseps
         correc = erlarg
-        ertest = dmax1(epsabs,epsrel*abs(reseps))
+        ertest = max(epsabs,epsrel*abs(reseps))
         if(abserr<=ertest) exit
       end if
 !
@@ -4314,7 +4314,7 @@ contains
 !       test on divergence
 !
     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
     if(0.1e-01_wp_>(result/area).or.(result/area)>0.1e+03_wp_ &
      .or.errsum>abs(area)) ier = 6
@@ -4416,7 +4416,7 @@ contains
     real(wp_), dimension(np), intent(in) :: apar
     real(wp_), intent(out) :: result,abserr,resabs,resasc
     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
     real(wp_), dimension(7) :: fv1,fv2
     integer :: j
@@ -4528,8 +4528,8 @@ contains
     resabs = resabs*hlgth
     abserr = abs((resk-resg)*hlgth)
     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_)
-    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = dmax1 &
+     min(0.1e+01_wp_,(0.2e+03_wp_*abserr/resasc)**1.5e+00_wp_)
+    if(resabs>uflow/(0.5e+02_wp_*epmach)) abserr = max &
      ((epmach*0.5e+02_wp_)*resabs,abserr)
   end subroutine dqk15imv