diff --git a/source/f90/boson.f90 b/source/f90/boson.f90
index c12f39c39dcafc66d5b67b47f0c4f993850a3b55..18f164bd17f4d1e1f1500da5d76cd827ee8cfdfb 100644
--- a/source/f90/boson.f90
+++ b/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
diff --git a/source/f90/doboson.f90 b/source/f90/doboson.f90
index 31a989b22b1d5056834efdc6c23e8cb568ac6f2f..7d53efa6aa91b4b8b3d66853b56a465670961376 100644
--- a/source/f90/doboson.f90
+++ b/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
diff --git a/source/f90/doeels.f90 b/source/f90/doeels.f90
index 56f28dafc49654e1d85607d4b051ad17d8362910..829415cc628348b473e08733354cf50f0ad683b0 100644
--- a/source/f90/doeels.f90
+++ b/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
 
diff --git a/source/f90/eels.f90 b/source/f90/eels.f90
index 99345eff336ef3809a613997ecc5513694456d0a..fb7a82d9d84d2f2cbc30da74c2cb55c9d6b589fc 100644
--- a/source/f90/eels.f90
+++ b/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
diff --git a/source/f90/fint3.f90 b/source/f90/fint3.f90
index 7135e6f590a0acc464a660e71f853402881ef047..35423688e1aa4e87119ff318a6fabf029e09da87 100644
--- a/source/f90/fint3.f90
+++ b/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
diff --git a/source/f90/quanc8.f90 b/source/f90/quanc8.f90
index 86ee4f36463a6f8150bb5f52b50528ca849879f5..a4a8dc6699480fc3ebf2b507dcb62ba21a3efd1e 100644
--- a/source/f90/quanc8.f90
+++ b/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)
diff --git a/source/f90/queels.f90 b/source/f90/queels.f90
index d2128df921e17e578004094b2218052d79c1cbe1..e40c6535a0821df894c2b34825ea8b004c66d808 100644
--- a/source/f90/queels.f90
+++ b/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 ***'
diff --git a/source/f90/rcffi.f90 b/source/f90/rcffi.f90
index aca42437a88babe690f32d5fdec7764e36e9f75b..e7db39b8e76614806ad635d4c88f3b8c514f4fe2 100644
--- a/source/f90/rcffi.f90
+++ b/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