refactor common / mulayr /

remove common / mulayr / by
- replacing the calculation of epsmac by call of function epsilon
- move calculation of argmin and argmax to the function surlos. They are used there only.

source/f90/doeels.f90

@@ -55,8 +55,8 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
 
   logical :: rational, user
   integer :: i, iw, neps, nofu, nout, nw, lmax
-  double precision :: a, acoef, aerr, alpha, argmin, argmax, b, bcoef, beta,  &
-      c1, c2, ccoef, cospsi, dlimf, dx, elleps, ener, epsmac, facru, f, f0,   &
+  double precision :: a, acoef, aerr, alpha, b, bcoef, beta,  &
+      c1, c2, ccoef, cospsi, dlimf, dx, elleps, ener, facru, f, f0,   &
       f1, fpic, fun, pi, prefac, psia, psii, qrat, rerr, ru, sinpsi, t,       &
       tanpsi, table, um, widt, wn, wpic, x, xmin, xmax, z, z1, z2
   double complex, allocatable :: eps(:)
@@ -66,7 +66,6 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
 
   common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
                    ru, um, dlimf, wn, user, rational
-  common / mulayr / argmin, argmax, epsmac
 
   data aerr / 0.0d0 /, rerr / 1.0d-06 /, f / 0.0d0 /, f1 / 0.0d0 /
 
@@ -86,12 +85,6 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
 ! *** epsmac + 1.0 = epsmac , cosh(argmin) = 1.0 , tanh(argmax) = 1.0
 
   pi = 4 * datan(1.0d0)
-  epsmac = 1.0d0
-  do while (1.0d0 + epsmac > 1.0d0)
-    epsmac = epsmac / 2
-  enddo
-  argmin = dsqrt(2 * epsmac)
-  argmax = 0.5d0 * dlog(2 / epsmac)
 
   dlimf = 0.0d0
   rational = .false.

source/f90/seteps.f90

@@ -31,12 +31,10 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
   double complex, intent(in out) :: eps(:)
   integer, intent(in) :: layers
 
-  double precision :: argmin, argmax, epsmac, x
+  double precision :: x
   double complex :: deno, nomi
   integer :: j, k, l, m
 
-  common / mulayr / argmin, argmax, epsmac
-
 !  write (*,*) 'seteps:'
 !  write (*,*) 'nos: ', size(nos)
 !  write (*,*) 'osc: ', size(osc)
@@ -78,8 +76,8 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode)
       if (osc(2, j) >= 0.0d0) then
         deno = 1.0d0 - deno
       endif
-      if (cdabs(deno) == 0.0d0) then
-        deno = epsmac
+      if (cdabs(deno) == 0.0d0) then ! replace 0 by machine epsilon
+        deno = epsilon(1.0d0) / 2
       endif
       eps(l) = eps(l) + osc(2, j) / deno
       enddo

source/f90/surlos.f90

@@ -28,7 +28,10 @@ double precision function surlos(dk, eps, thick, layers, nper)
   double precision :: argmin, argmax, cn, cnm1, epsmac, sn, snm1, tn
   double complex :: a, b, csi, pnm2, pnm1, pn, pp, qnm2, qnm1, qn, qp, z
 
-  common / mulayr / argmin, argmax, epsmac
+  epsmac = epsilon(1.0d0)
+  argmin = dsqrt(2 * epsmac)
+  epsmac = epsmac / 2
+  argmax = 0.5d0 * dlog(2 / epsmac)
 
   zero(z) = (dble(z) == 0.0d0) .and. (dimag(z) == 0.0d0)