source/f90/boson.f90

@@ -50,16 +50,26 @@ program boson
   np = 0
   do
     read(12, *, IOSTAT = ioStatus) wmax, y
-    if (ioStatus /= 0) exit
-    if (wmax < 0.0d0) cycle
+    if (ioStatus /= 0) then
+      exit
+    endif
+    if (wmax < 0.0d0) then
+      cycle
+    endif
     np = np + 1
     p2(np) = y
-    if (np == 1) wmin = wmax
-    if (np < nmax) cycle
+    if (np == 1) then
+      wmin = wmax
+    endif
+    if (np < nmax) then
+      cycle
+    endif
   enddo
   close(unit = 12)
 
-  if (np <= 0) stop '*** no input values for pcl ***'
+  if (np <= 0) then
+    stop '*** no input values for pcl ***'
+  endif
 
   write(*,*) comment(2)
   write(*,'(a, i6, a, g15.7, a, g15.7)') ' read', np, ' values of pcl from', wmin, ' to', wmax

source/f90/doboson.f90

@@ -76,7 +76,9 @@ subroutine doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p, debug,
   do m = 2, mmax
     n = 2 * n
     wmax = n * dwn
-    if (wmax >= test) exit
+    if (wmax >= test) then
+      exit
+    endif
   enddo
   if (wmax < test) then
     m = mmax
@@ -226,9 +228,13 @@ subroutine doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p, debug,
     endif
     do i = imin, n
       j = n - i
-      if (mod(j, istep) > 0) cycle
+      if (mod(j, istep) > 0) then
+        cycle
+      endif
       x = -j * dwn
-      if (x > emax) exit
+      if (x > emax) then
+        exit
+      endif
       nout = nout + 1
       xout(nout)= x
       yout(nout)= p2(j + 1)
@@ -242,9 +248,13 @@ subroutine doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p, debug,
       endif
       if (imax >= 2) then
         do i = 2, imax
-          if (mod(i - 1, istep) > 0) cycle
+          if (mod(i - 1, istep) > 0) then
+            cycle
+          endif
           x = (i - 1) * dwn
-          if (x < emin) cycle
+          if (x < emin) then
+            cycle
+          endif
           nout = nout + 1
           xout(nout) = x
           yout(nout) = p1(i)
@@ -272,7 +282,9 @@ subroutine doboson(t, width, gauss, asym, emin, emax, wmin, wmax, np, p, debug,
     fp1 = fp
     fm = p2(i)
     fp = p1(i)
-    if (i == 2) cycle
+    if (i == 2) then
+      cycle
+    endif
     picm = .false.
     picp = .false.
     wn = (i - 1) * dwn

source/f90/doeels.f90

@@ -91,7 +91,9 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
   user = layers == 0
   if (user) then
 
-    if (layers == 1) rational = .true.
+    if (layers == 1) then
+      rational = .true.
+    endif
     if (contrl == 'image') then
 ! *** image-charge screening factor
       if (layers == 1 .and. neps == 2) then
@@ -150,7 +152,9 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
       c1 = c1 * c2
       xmin = wmin / (2 * ener * psia)
       xmax = wmax / (2 * ener * psia)
-      if (xmin <= 0.0d0) xmin = 0.0d0
+      if (xmin <= 0.0d0) then
+        xmin = 0.0d0
+      endif
       dx = dmax1(0.02d0, (xmax - xmin) / nt)
       z1 = 0.0d0
       z2 = 0.0d0
@@ -159,9 +163,13 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
         call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
         table(i) = f
         f = f * (1.0d0 + alpha * x)**2
-        if (dabs(c2 * f - c1) < c2 * rerr) cycle
+        if (dabs(c2 * f - c1) < c2 * rerr) then
+          cycle
+        endif
         z = (1.0d0 - f) / (c2 * f - c1)
-        if (z <= 0.0d0) cycle
+        if (z <= 0.0d0) then
+          cycle
+        endif
         z1 = z1 + x * z * (x**2 - z)
         z2 = z2 + (x * z)**2
       enddo
@@ -196,18 +204,7 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
     wn = wmin + (iw - 1) * dw
 !    if (debug) write (*,*) 'wn: ', wn
     if (wn >= 0.0d0) then
-      if (wn /= 0.0d0) then
-        if (.not. user) call seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
-
-        x = wn / (2 * ener * psia)
-        if (rational) then
-          f = qrat(x, alpha, beta, c1, c2) * dimag(-2 / (1.0d0 + eps(1)))
-        else
-          call queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
-        endif
-        f = prefac * f / wn
-      endif ! wn /= 0.0d0
-
+q
       wn_array(iw) = wn
       f_array(iw)  = f
 

source/f90/eels.f90

@@ -88,8 +88,12 @@ program eels
       read(11, *) epsinf(l), nos(l)
       write(6, 103)
       if (nos(l) <= 0) then
-        if (l <= layers) write(6, 104) l, name(l), thick(l), epsinf(l)
-        if (l  > layers) write(6, 105) epsinf(l)
+        if (l <= layers) then
+          write(6, 104) l, name(l), thick(l), epsinf(l)
+        endif
+        if (l  > layers) then
+          write(6, 105) epsinf(l)
+        endif
       else
         do j = 1, nos(l)
           if (.not. allocated(osc)) then

source/f90/fint3.f90

@@ -37,14 +37,20 @@ double precision function fint3(u, eps, thick, layers, nper, eps_size)
   endif
   rac = dsign(1.0d0, acoef) * cospsi * dsqrt((1.0d0 - elleps) * acoef * (um - u) * (um + u))
   arg = (bcoef - rac) / (u * acoef)
-  if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg)
+  if (dabs(arg) > 1.0d0) then
+    arg = dsign(1.0d0, arg)
+  endif
   phi2 = dacos(arg)
   fint3 = phint(phi2, tanpsi, u)
   arg = (bcoef + rac) / (u * acoef)
-  if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg)
+  if (dabs(arg) > 1.0d0) then
+    arg = dsign(1.0d0, arg)
+  endif
   phi1 = dacos(arg)
   fint3 = fint3 - phint(phi1, tanpsi, u)
-  if (rational) return
+  if (rational) then
+    return
+  endif
   if (user) then
     fint3 = fint3 * usurlo(ru * u, wn)
   else

source/f90/quanc8.f90

@@ -88,7 +88,9 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, t
   errest = 0.0d0
   area   = 0.0d0
   nofun = 0
-  if (a == b) return
+  if (a == b) then
+    return
+  endif
 
 ! ***   stage 2 ***   initialization for first interval
 
@@ -193,7 +195,9 @@ subroutine quanc8(fun, a, b, abserr, relerr, result, errest, nofun, flag, eps, t
       enddo
       nim = nim + 1
 
-      if (lev <= 0) exit
+      if (lev <= 0) then
+        exit
+      endif
 
 ! assemble elements required for the next interval.
       qprev = qright(lev)

source/f90/queels.f90

@@ -80,7 +80,9 @@ subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
     flag(3) = 0.0d0
   endif
   do ie = 1, 3
-    if (flag(ie) == 0.0d0) cycle
+    if (flag(ie) == 0.0d0) then
+      cycle
+    endif
     write(*,*) ' +++ flag(', ie, ') =', flag(ie), ', error =', error(ie), ' +++'
     if (flag(ie) - aint(flag(ie)) > 0.5d-02) then
       stop '*** execution aborted ***'

source/f90/rcffi.f90

@@ -104,7 +104,9 @@ subroutine rcffi(ar, ai, msign, h)
     ai(i + 1) = ai(i) - ti
     ai(i) = ai(i) + ti
   enddo
-  if (m == 1) return
+  if (m == 1) then
+    return
+  endif
   c = 0.0d0
   s = dble(isign(1, msign))
   le = 2

source/f90/seteps.f90

@@ -52,8 +52,8 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
           yTO = osc(2,k)
           wLO = osc(3,k)
           yLO = osc(4,k)
-          nomi = nomi * dcmplx(wLO**2 - wn2, yLO * wn)
-          deno = deno * dcmplx(wTO**2 - wn2, yTO * wn)
+          nomi = nomi * dcmplx(wLO**2 - wn2, -yLO * wn)
+          deno = deno * dcmplx(wTO**2 - wn2, -yTO * wn)
         case (Drude)! Lambin 1990, equ. 8
           wp = osc(1,k)
           y  = osc(2,k)
@@ -62,7 +62,7 @@ subroutine seteps(neps, nos, oscType, osc, epsinf, wn, eps, layers)
           wp = osc(1,k)
           y0 = osc(2,k)
           yp = osc(3,k)
-          addDrude = addDrude - dcmplx(wp**2, (yp - y0) * wn) / dcmplx(wn2, wn * y0)
+          addDrude = addDrude - dcmplx(wp**2, (y0 - yp) * wn) / dcmplx(wn2, wn * y0)
         case default
           write (*,'(A)') '*** Error in oscillator type! ***'
       end select