common param replaced by module

source/f90/Makefile

@@ -4,7 +4,8 @@
 FC = gfortran 
 
 # fortran compiler options
-FFLAGS = -g -gdwarf-2 -fbounds-check -fcheck=all
+FFLAGS = -g -gdwarf-2 -fbounds-check -fcheck=all -ffpe-trap=invalid -O0 -Wall
+# -Wextra -Werror -Wpedantic
 # gfortran version 4.8 does not know -fdiagnostics-color
 # safeguard for major version >= 5
 GFORTAN_VERSION_GE_5 := $(shell echo `gcc -dumpversion | cut -f1 -d. ` \>= 5 | bc)
@@ -30,14 +31,15 @@ bosonsubs = doboson.o respon.o sicot.o sintr.o rcffi.o o.o
 boson: boson.f90 change_working_dir.o $(bosonsubs) $(bosonmods)
 	$(FC) $(FFLAGS) -o bosonf90 boson.f90 change_working_dir.o $(bosonsubs) 
 
-eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod
+eelsmods = quanc8_mod.mod queels_mod.mod seteps_mod.mod param_mod.mod
 doeels.o: doeels.f90 $(eelsmods)
 
 queels_mod.mod: quanc8_mod.mod fint1_mod.mod fint2_mod.mod fint3_mod.mod
 quanc8_mod.mod: fun_mod.mod
-fint1_mod.mod: surlos_mod.mod
-fint2_mod.mod: surlos_mod.mod
-fint3_mod.mod: surlos_mod.mod
+fint1_mod.mod: surlos_mod.mod param_mod.mod
+fint2_mod.mod: surlos_mod.mod param_mod.mod
+fint3_mod.mod: surlos_mod.mod param_mod.mod
+fun_mod.mod: param_mod.mod
 
 getoptsubs = getopt.o text.o system.o constants.o date_and_time.o kinds.o dummy_variables.o
 
@@ -64,7 +66,7 @@ sufr_kinds.mod:
 sufr_dummy.mod:
 	$(FC) $(FFLAGS) -c -o dummy_variables.o ../getopt-libs/libsufr-0.7.7/src/dummy_variables.f90
 
-eelssubs = doeels.o usurlo.o quanc8.o fun.o queels.o fint1.o fint2.o fint3.o surlos.o seteps.o phint.o qrat.o
+eelssubs = doeels.o usurlo.o quanc8.o fun.o queels.o fint1.o fint2.o fint3.o surlos.o seteps.o phint.o qrat.o param.o
 eels: eels.f90 change_working_dir.o $(eelssubs) $(eelsmods)
 	$(FC) $(FFLAGS) -o eelsf90 eels.f90 change_working_dir.o $(eelssubs)
 

source/f90/doeels.f90

@@ -38,6 +38,7 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
   use quanc8_mod
   use queels_mod
   use seteps_mod
+  use param_mod
 
   implicit none
 
@@ -53,20 +54,16 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
   double precision, intent(in out) :: wn_array(wn_array_size), f_array(wn_array_size)
   logical, intent(in), optional :: debug
 
-  logical :: rational, user
   integer :: i, iw, neps, nofu, nout, nw, lmax
-  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 precision :: a, aerr, alpha, b, beta,              &
+      c1, c2, dx, ener, facru, f, f0,                       &
+      f1, fpic, fun, pi, prefac, psia, psii, qrat, rerr, t, &
+      width, wpic, x, xmin, xmax, z, z1, z2
+  double precision :: table(nt)
   double complex, allocatable :: eps(:)
-  dimension table(nt)
 
   external fun, qrat
 
-  common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
-                   ru, um, dlimf, wn, user, rational
-
   data aerr / 0.0d0 /, rerr / 1.0d-06 /, f / 0.0d0 /, f1 / 0.0d0 /
 
   if (debug) then
@@ -218,8 +215,8 @@ subroutine doeels (e0, theta, phia, phib, wmin, wmax, dw, comment, comment_size,
           t = b / a
           wpic = wn - t * dw
           fpic = f + 0.5d0 * b * t
-          widt = dsqrt(8 * fpic / a) * dw
-          if (debug) write(*, 120) wpic, fpic, widt
+          width = dsqrt(8 * fpic / a) * dw
+          if (debug) write(*, 120) wpic, fpic, width
         endif ! a > 4 * rerr * f1
       endif ! iw >= 3 ...
     endif ! wn >= 0.0d0

source/f90/fint1.f90

@@ -15,6 +15,7 @@ double precision function fint1(u, eps, thick, layers, nper, eps_size)
 ! ******************************************************************
 
   use surlos_mod
+  use param_mod
 
   implicit none
 
@@ -23,13 +24,8 @@ double precision function fint1(u, eps, thick, layers, nper, eps_size)
   double precision, intent(in) :: thick(eps_size)
   integer, intent(in) :: layers, nper, eps_size
 
-  logical :: rational, user
-  double precision :: acoef, bcoef, ccoef, cospsi, den, dif, dlimf, e, elleps
-  double precision :: pi, rom, rop, ru, sinpsi, sum, um, usurlo, t
-  double precision :: tanpsi, wn, u2
-
-  common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
-      ru, um, dlimf, wn, user, rational
+  double precision :: den, dif, e, pi, rom, rop, sum, t, u2
+  double precision :: usurlo
 
   data pi / 3.141592653589793238d0 /
 

source/f90/fint2.f90

@@ -15,6 +15,7 @@ double precision function fint2(u, eps, thick, layers, nper, eps_size)
 ! ******************************************************************
 
   use surlos_mod
+  use param_mod
 
   implicit none
 
@@ -23,12 +24,8 @@ double precision function fint2(u, eps, thick, layers, nper, eps_size)
   double precision, intent(in) :: thick(eps_size)
   integer, intent(in) :: layers, nper, eps_size
 
-  logical :: rational, user
-  double precision :: a, arg, b, b2, c, ccoef, cospsi, dlimf, elleps, phi
-  double precision :: phint, ru, sinpsi, um, usurlo, t, tanpsi, wn, x
-
-  common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi, &
-      ru, um, dlimf, wn, user, rational
+  double precision :: arg, b2, c, phi, t, x
+  double precision :: phint, usurlo
 
 !  write (*,*) 'fint2:'
 !  write (*,*) 'thick: ', size(thick)
@@ -38,13 +35,13 @@ double precision function fint2(u, eps, thick, layers, nper, eps_size)
     fint2 = 0.0d0
     return
   endif
-  b2 = b**2
-  c = (1.0d0 - elleps) * (cospsi * u)**2 + (b - ccoef) * (b + ccoef)
-  if (dabs(a * c) > 1.0d-03 * b2) then
-    x = (b - dsqrt(b2 - a * c)) / a
+  b2 = bcoef**2
+  c = (1.0d0 - elleps) * (cospsi * u)**2 + (bcoef - ccoef) * (bcoef + ccoef)
+  if (dabs(acoef * c) > 1.0d-03 * b2) then
+    x = (bcoef - dsqrt(b2 - acoef * c)) / acoef
   else
-    x = a * c / b2
-    x = 0.5d0 * c * (1.d0 + 0.25d0 * x * (1.d0 + 0.5d0 * x * (1.d0 + 0.625d0 * x))) / b
+    x = acoef * c / b2
+    x = 0.5d0 * c * (1.d0 + 0.25d0 * x * (1.d0 + 0.5d0 * x * (1.d0 + 0.625d0 * x))) / bcoef
   endif
   arg = x / u
   if (dabs(arg) > 1.0d0) then

source/f90/fint3.f90

@@ -15,6 +15,7 @@ double precision function fint3(u, eps, thick, layers, nper, eps_size)
 ! ******************************************************************
 
   use surlos_mod
+  use param_mod
 
   implicit none
 
@@ -23,12 +24,8 @@ double precision function fint3(u, eps, thick, layers, nper, eps_size)
   double precision, intent(in) :: thick(eps_size)
   integer, intent(in) :: layers, nper, eps_size
 
-  logical :: rational, user
-  double precision :: a, arg, b, ccoef, cospsi, dlimf, elleps, phi1, phi2
-  double precision :: phint, sinpsi, rac, ru, um, usurlo, t, tanpsi, wn
-
-  common / param / a, b, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
-      ru, um, dlimf, wn, user, rational
+  double precision :: arg, phi1, phi2, rac, t
+  double precision :: phint, usurlo
 
 !  write (*,*) 'fint3:'
 !  write (*,*) 'thick: ', size(thick)
@@ -38,12 +35,12 @@ double precision function fint3(u, eps, thick, layers, nper, eps_size)
     fint3 = 0.0d0
     return
   endif
-  rac = dsign(1.0d0, a) * cospsi * dsqrt((1.0d0 - elleps) * a * (um - u) * (um + u))
-  arg = (b - rac) / (u * a)
+  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)
   phi2 = dacos(arg)
   fint3 = phint(phi2, tanpsi, u)
-  arg = (b + rac) / (u * a)
+  arg = (bcoef + rac) / (u * acoef)
   if (dabs(arg) > 1.0d0) arg = dsign(1.0d0, arg)
   phi1 = dacos(arg)
   fint3 = fint3 - phint(phi1, tanpsi, u)

source/f90/fun.f90

@@ -7,16 +7,13 @@ double precision function fun(phi)
 ! *                                                                *
 ! ******************************************************************
 
+  use param_mod
+
   implicit none
 
   double precision, intent(in) :: phi
 
-  logical :: user, rational
-  double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps, ru
-  double precision :: sinphi, sinpsi, tanpsi, um, wn
-
-  common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
-                   ru, um, dlimf, wn, user, rational
+  double precision :: sinphi
 
   sinphi = dsin(phi)
   fun = dsqrt((1.0d0 - elleps + elleps * sinphi**2) *   &

source/f90/param.f90

+module param_mod
+  double precision :: acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi
+  double precision :: ru, um, dlimf, wn
+  logical :: user, rational
+end module param_mod

source/f90/queels.f90

@@ -23,6 +23,7 @@ subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
   use fint1_mod
   use fint2_mod
   use fint3_mod
+  use param_mod
 
   implicit none
 
@@ -35,16 +36,11 @@ subroutine queels(x, f, aerr, rerr, facru, eps, thick, layers, nper)
   double precision, intent(in) :: thick(:)
   integer, intent(in) :: layers, nper
 
-  logical :: rational, user
-  double precision :: acoef, bcoef, ccoef, cospsi, dlimf, elleps
-  double precision :: error, flag, ru, sinpsi
-  double precision :: u1, u2, um, ut, tanpsi, wn, y
+  double precision :: error, flag
+  double precision :: u1, u2, ut, y
   integer :: ie, nofu
   dimension error(3), flag(3)
 
-  common / param / acoef, bcoef, ccoef, elleps, cospsi, sinpsi, tanpsi,  &
-      ru, um, dlimf, wn, user, rational
-
 !  write (*,*) 'queels:'
 !  write (*,*) 'eps: ', size(eps)
 !  write (*,*) 'thick: ', size(thick)

source/f90/surlos.f90

@@ -32,7 +32,7 @@ double precision function surlos(dk, eps, thick, layers, nper)
   integer :: lstart, n
   double precision, allocatable :: arg(:)
   double precision :: argmin, argmax, cn, cnm1, epsmac, sn, snm1, tn
-  double complex :: a, b, csi, pnm2, pnm1, pn, pp, qnm2, qnm1, qn, qp, z
+  double complex :: a, b, csi, pnm2, pnm1, pn, pp, qnm2, qnm1, qn, qp
 
   epsmac = epsilon(1.0d0)
   argmin = dsqrt(2 * epsmac)