Fix more instances of double values, remove instances of dabs, dsqrt,…

… and dexp, improve m_precision_select, incorporate changes from patch file

src/common/m_constants.fpp

@@ -8,11 +8,11 @@ module m_constants
 
     character, parameter :: dflt_char = ' ' !< Default string value
 
-    real(wp), parameter :: dflt_real = -1d6                !< Default real value
-    real(wp), parameter :: sgm_eps = 1d-16               !< Segmentation tolerance
-    real(wp), parameter :: small_alf = 1d-11                !< Small alf tolerance
+    real(wp), parameter :: dflt_real = (-1_wp*(10._wp**6))                !< Default real value
+    real(wp), parameter :: sgm_eps = (1_wp*(10._wp**-16))               !< Segmentation tolerance
+    real(wp), parameter :: small_alf = (1_wp*(10._wp**-11))                !< Small alf tolerance
     real(wp), parameter :: pi = 3.141592653589793_wp !< Pi
-    real(wp), parameter :: verysmall = 1.d-12              !< Very small number
+    real(wp), parameter :: verysmall = (1._wp*(10._wp**-12))              !< Very small number
 
     integer, parameter :: num_stcls_min = 5    !< Minimum # of stencils
     integer, parameter :: path_len = 400  !< Maximum path length

src/common/m_eigen_solver.f90

@@ -163,8 +163,8 @@ subroutine cbal(nm, nl, ar, ai, low, igh, scale)
 
             do 200 j = k, l
                 if (j == i) go to 200
-                c = c + dabs(ar(j, i)) + dabs(ai(j, i))
-                r = r + dabs(ar(i, j)) + dabs(ai(i, j))
+                c = c + abs(ar(j, i)) + abs(ai(j, i))
+                r = r + abs(ar(i, j)) + abs(ai(i, j))
 200         end do
 !     .......... guard against zero c or r due to underflow ..........
             if (c == 0.0_wp .or. r == 0.0_wp) go to 270
@@ -243,7 +243,7 @@ subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
             scale = 0.0_wp
 !     .......... scale column (algol tol then not needed) ..........
             do 90 i = ml, igh
-                scale = scale + dabs(ar(i, ml - 1)) + dabs(ai(i, ml - 1))
+                scale = scale + abs(ar(i, ml - 1)) + abs(ai(i, ml - 1))
 90          end do
             if (scale == 0._wp) go to 180
             mp = ml + igh
@@ -255,7 +255,7 @@ subroutine corth(nm, nl, low, igh, ar, ai, ortr, orti)
                 h = h + ortr(i)*ortr(i) + orti(i)*orti(i)
 100         end do
 !
-            g = dsqrt(h)
+            g = sqrt(h)
             call pythag(ortr(ml), orti(ml), f)
             if (f == 0._wp) go to 103
             h = h + f*g
@@ -375,8 +375,8 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 !     .......... for i=igh-1 step -1 until low+1 do -- ..........
 105     do 140 ii = 1, iend
             i = igh - ii
-            if (dabs(ortr(i)) == 0._wp .and. dabs(orti(i)) == 0._wp) go to 140
-            if (dabs(hr(i, i - 1)) == 0._wp .and. dabs(hi(i, i - 1)) == 0._wp) go to 140
+            if (abs(ortr(i)) == 0._wp .and. abs(orti(i)) == 0._wp) go to 140
+            if (abs(hr(i, i - 1)) == 0._wp .and. abs(hi(i, i - 1)) == 0._wp) go to 140
 !     .......... norm below is negative of h formed in corth ..........
             norm = hr(i, i - 1)*ortr(i) + hi(i, i - 1)*orti(i)
             ip1 = i + 1
@@ -411,7 +411,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 !
         do 170 i = l, igh
             ll = min0(i + 1, igh)
-            if (dabs(hi(i, i - 1)) == 0._wp) go to 170
+            if (abs(hi(i, i - 1)) == 0._wp) go to 170
             call pythag(hr(i, i - 1), hi(i, i - 1), norm)
             yr = hr(i, i - 1)/norm
             yi = hi(i, i - 1)/norm
@@ -456,9 +456,9 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 240     do 260 ll = low, en
             l = en + low - ll
             if (l == low) go to 300
-            tst1 = dabs(hr(l - 1, l - 1)) + dabs(hi(l - 1, l - 1)) &
-                   + dabs(hr(l, l)) + dabs(hi(l, l))
-            tst2 = tst1 + dabs(hr(l, l - 1))
+            tst1 = abs(hr(l - 1, l - 1)) + abs(hi(l - 1, l - 1)) &
+                   + abs(hr(l, l)) + abs(hi(l, l))
+            tst2 = tst1 + abs(hr(l, l - 1))
             if (tst2 == tst1) go to 300
 260     end do
 !     .......... form shift ..........
@@ -481,7 +481,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
         si = si - xxi
         go to 340
 !     .......... form exceptional shift ..........
-320     sr = dabs(hr(en, enm1)) + dabs(hr(enm1, en - 2))
+320     sr = abs(hr(en, enm1)) + abs(hr(enm1, en - 2))
         si = 0.0_wp
 !
 340     do 360 i = low, en
@@ -523,7 +523,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 500     end do
 !
         si = hi(en, en)
-        if (dabs(si) == 0._wp) go to 540
+        if (abs(si) == 0._wp) go to 540
         call pythag(hr(en, en), si, norm)
         sr = hr(en, en)/norm
         si = si/norm
@@ -568,7 +568,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 590         end do
 600     end do
 !
-        if (dabs(si) == 0._wp) go to 240
+        if (abs(si) == 0._wp) go to 240
 !
         do 630 i = 1, en
             yr = hr(i, en)
@@ -598,7 +598,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 !
         do i = 1, nl
             do j = i, nl
-                tr = dabs(hr(i, j)) + dabs(hi(i, j))
+                tr = abs(hr(i, j)) + abs(hi(i, j))
                 if (tr > norm) norm = tr
             end do
         end do
@@ -635,7 +635,7 @@ subroutine comqr2(nm, nl, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr)
 765             continue
                 call cdiv(zzr, zzi, yr, yi, hr(i, en), hi(i, en))
 !     .......... overflow control ..........
-                tr = dabs(hr(i, en)) + dabs(hi(i, en))
+                tr = abs(hr(i, en)) + abs(hi(i, en))
                 if (tr == 0.0_wp) go to 780
                 tst1 = tr
                 tst2 = tst1 + 1.0_wp/tst1
@@ -709,12 +709,12 @@ end subroutine comqr2
         !!           transformed in their first ml columns
     subroutine cbabk2(nm, nl, low, igh, scale, ml, zr, zi)
         integer, intent(in) :: nm, nl, low, igh
-        double precision, intent(in) :: scale(nl)
+        real(wp) :: scale(nl)
         integer, intent(in) :: ml
-        double precision, intent(inout) :: zr(nm, ml), zi(nm, ml)
+        real(wp) :: zr(nm, ml), zi(nm, ml)
 
         integer :: i, j, k, ii
-        double precision :: s
+        real(wp) :: s
 
         if (ml == 0) go to 200
         if (igh == low) go to 120
@@ -757,14 +757,14 @@ subroutine csroot(xr, xi, yr, yi)
         real(wp), intent(in) :: xr, xi
         real(wp), intent(out) :: yr, yi
 !
-!     (yr,yi) = complex dsqrt(xr,xi)
+!     (yr,yi) = complex sqrt(xr,xi)
 !     branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi)
 !
         real(wp) :: s, tr, ti, c
         tr = xr
         ti = xi
         call pythag(tr, ti, c)
-        s = dsqrt(0.5_wp*(c + dabs(tr)))
+        s = sqrt(0.5_wp*(c + abs(tr)))
         if (tr >= 0.0_wp) yr = s
         if (ti < 0.0_wp) s = -s
         if (tr <= 0.0_wp) yi = s
@@ -786,7 +786,7 @@ subroutine cdiv(ar, ai, br, bi, cr, ci)
         ! cr = (ar*br + ai*bi) / (br**2._wp + bi**2._wp)
         ! ci = (ai*br - ar*bi) / (br**2._wp + bi**2._wp)
 
-        s = dabs(br) + dabs(bi)
+        s = abs(br) + abs(bi)
         ars = ar/s
         ais = ai/s
         brs = br/s
@@ -801,12 +801,12 @@ subroutine pythag(a, b, c)
         real(wp), intent(in) :: a, b
         real(wp), intent(out) :: c
 !
-!     finds dsqrt(a**2+b**2) without overflow or destructive underflow
+!     finds sqrt(a**2+b**2) without overflow or destructive underflow
 !
         real(wp) :: p, r, s, t, u
-        p = dmax1(dabs(a), dabs(b))
+        p = dmax1(abs(a), abs(b))
         if (p == 0.0_wp) go to 20
-        r = (dmin1(dabs(a), dabs(b))/p)**2
+        r = (dmin1(abs(a), abs(b))/p)**2
 10      continue
         t = 4.0_wp + r
         if (t == 4.0_wp) go to 20

src/common/m_helper.fpp

@@ -68,7 +68,7 @@ contains
         real(wp) :: nR3
 
         nR3 = dot_product(weights, nRtmp**3._wp)
-        ntmp = DSQRT((4._wp*pi/3._wp)*nR3/vftmp)
+        ntmp = sqrt((4._wp*pi/3._wp)*nR3/vftmp)
         !ntmp = (3._wp/(4._wp*pi))*0.00001
 
         !print *, "nbub", ntmp
@@ -153,8 +153,8 @@ contains
         if (thermal == 2) gamma_m = 1._wp
 
         temp = 293.15_wp
-        D_m = 0.242d-4
-        uu = DSQRT(pl0/rhol0)
+        D_m = (0.242_wp*(10._wp**-4))
+        uu = sqrt(pl0/rhol0)
 
         omega_ref = 3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/Web
 
@@ -163,10 +163,10 @@ contains
         R_n = Ru/M_n
         R_v = Ru/M_v
         ! phi_vn & phi_nv (phi_nn = phi_vv = 1)
-        phi_vn = (1._wp + DSQRT(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
-                 /(DSQRT(8._wp)*DSQRT(1._wp + M_v/M_n))
-        phi_nv = (1._wp + DSQRT(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
-                 /(DSQRT(8._wp)*DSQRT(1._wp + M_n/M_v))
+        phi_vn = (1._wp + sqrt(mu_v/mu_n)*(M_n/M_v)**(0.25_wp))**2 &
+                 /(sqrt(8._wp)*sqrt(1._wp + M_v/M_n))
+        phi_nv = (1._wp + sqrt(mu_n/mu_v)*(M_v/M_n)**(0.25_wp))**2 &
+                 /(sqrt(8._wp)*sqrt(1._wp + M_n/M_v))
         ! internal bubble pressure
         pb0 = pl0 + 2._wp*ss/(R0ref*R0)
 
@@ -208,7 +208,7 @@ contains
         !end if
 
         ! natural frequencies
-        omegaN = DSQRT(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
+        omegaN = sqrt(3._wp*k_poly*Ca + 2._wp*(3._wp*k_poly - 1._wp)/(Web*R0))/R0
         do ir = 1, Nb
             call s_transcoeff(omegaN(ir)*R0(ir), Pe_T(ir)*R0(ir), &
                               Re_trans_T(ir), Im_trans_T(ir))
@@ -260,43 +260,43 @@ contains
         real(wp), dimension(nb) :: phi
 
         ! nondiml. min. & max. initial radii for numerical quadrature
-        !sd   = 0.05D0
-        !R0mn = 0.75D0
-        !R0mx = 1.3D0
+        !sd   = (0.05_wp * (10._wp ** 0))
+        !R0mn = (0.75_wp * (10._wp ** 0))
+        !R0mx = (1.3_wp * (10._wp ** 0))
 
-        !sd   = 0.3D0
-        !R0mn = 0.3D0
-        !R0mx = 6.D0
+        !sd   = (0.3_wp * (10._wp ** 0))
+        !R0mn = (0.3_wp * (10._wp ** 0))
+        !R0mx = (6._wp * (10._wp ** 0))
 
-        !sd   = 0.7D0
-        !R0mn = 0.12D0
-        !R0mx = 150.D0
+        !sd   = (0.7_wp * (10._wp ** 0))
+        !R0mn = (0.12_wp * (10._wp ** 0))
+        !R0mx = (150._wp * (10._wp ** 0))
 
         sd = poly_sigma
-        R0mn = 0.8_wp*DEXP(-2.8_wp*sd)
-        R0mx = 0.2_wp*DEXP(9.5_wp*sd) + 1._wp
+        R0mn = 0.8_wp*exp(-2.8_wp*sd)
+        R0mx = 0.2_wp*exp(9.5_wp*sd) + 1._wp
 
         ! phi = ln( R0 ) & return R0
         do ir = 1, nb
-            phi(ir) = DLOG(R0mn) &
-                      + dble(ir - 1)*DLOG(R0mx/R0mn)/dble(nb - 1)
-            R0(ir) = DEXP(phi(ir))
+            phi(ir) = log(R0mn) &
+                      + dble(ir - 1)*log(R0mx/R0mn)/dble(nb - 1)
+            R0(ir) = exp(phi(ir))
         end do
         dphi = phi(2) - phi(1)
 
         ! weights for quadrature using Simpson's rule
         do ir = 2, nb - 1
             ! Gaussian
-            tmp = DEXP(-0.5_wp*(phi(ir)/sd)**2)/DSQRT(2._wp*pi)/sd
+            tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
             if (mod(ir, 2) == 0) then
                 weight(ir) = tmp*4._wp*dphi/3._wp
             else
                 weight(ir) = tmp*2._wp*dphi/3._wp
             end if
         end do
-        tmp = DEXP(-0.5_wp*(phi(1)/sd)**2)/DSQRT(2._wp*pi)/sd
+        tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
         weight(1) = tmp*dphi/3._wp
-        tmp = DEXP(-0.5_wp*(phi(nb)/sd)**2)/DSQRT(2._wp*pi)/sd
+        tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
         weight(nb) = tmp*dphi/3._wp
     end subroutine s_simpson
 

src/common/m_helper_basic.f90

@@ -22,7 +22,7 @@ module m_helper_basic
     !> This procedure checks if two floating point numbers of wp are within tolerance.
     !! @param a First number.
     !! @param b Second number.
-    !! @param tol_input Relative error (default = 1d-6).
+    !! @param tol_input Relative error (default = (1_wp * (10._wp ** -6))).
     !! @return Result of the comparison.
     logical function f_approx_equal(a, b, tol_input) result(res)
         !$acc routine seq
@@ -35,7 +35,7 @@ logical function f_approx_equal(a, b, tol_input) result(res)
         if (present(tol_input)) then
             tol = tol_input
         else
-            tol = 1d-6
+            tol = (1_wp*(10._wp**-6))
         end if
 
         if (a == b) then