diff --git a/source/f90/seteps.f90 b/source/f90/seteps.f90 index 743ff75f721b926faf9371e65dc40147b984df7f..416d0cc99d0782d2b69410a30ef7816ece7da4e1 100644 --- a/source/f90/seteps.f90 +++ b/source/f90/seteps.f90 @@ -70,20 +70,20 @@ subroutine seteps(neps, nos, osc, epsinf, wn, name, eps, layers, mode) ! The original version had this additional loop. It seems, it has been removed ! because all cases of nos(l) > 1 are treated in case 1 above do k = 1, nos(l) - j = j + 1 - if (osc(1, j) > 1d-3) then ! prevent division by zero - x = wn / osc(1, j) - else - x = wn * 1d3 - endif - deno = x * dcmplx(x, osc(3, j)) - if (osc(2, j) >= 0.0d0) then - deno = 1.0d0 - deno - endif - 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 + j = j + 1 + if (osc(1, j) > 1d-3) then ! prevent division by zero + x = wn / osc(1, j) + else + x = wn * 1d3 + endif + deno = x * dcmplx(x, osc(3, j)) + if (osc(2, j) >= 0.0d0) then + deno = 1.0d0 - deno + endif + 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 endif enddo