diff --git a/src/calcei_mod.f90 b/src/calcei_mod.f90
index 4413918..e57f502 100644
--- a/src/calcei_mod.f90
+++ b/src/calcei_mod.f90
@@ -1,6 +1,8 @@
 module calcei_mod
 
- contains
+  implicit none
+
+contains
 
 ! ======================================================================
 ! nist guide to available math software.
diff --git a/src/const_and_precisions.f90 b/src/const_and_precisions.f90
index 048ae6c..d9d0fdb 100755
--- a/src/const_and_precisions.f90
+++ b/src/const_and_precisions.f90
@@ -34,6 +34,7 @@
   REAL(wp_), PARAMETER :: zero    = 0.0_wp_
   REAL(wp_), PARAMETER :: one     = 1.0_wp_
   real(wp_), parameter :: pi      = 3.141592653589793_wp_  ! 3.141592653589793238462643383280
+  real(wp_), parameter :: pihalf  = 1.57079632679489661923_wp_
   REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
 ! REAL(wp_), PARAMETER :: sqrt_2  = 1.414213562373095_wp_
   REAL(wp_), PARAMETER :: degree  = pi/180.0_wp_
diff --git a/src/dispersion.f90 b/src/dispersion.f90
index 1c9876d..34c8218 100644
--- a/src/dispersion.f90
+++ b/src/dispersion.f90
@@ -1,7 +1,7 @@
 module dispersion
 !
   use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi
-  use calcei_mod
+  use calcei_mod, only : calcei3
   implicit none
 ! local constants
   integer, parameter :: npts=500
diff --git a/src/dqagmv.f b/src/dqagmv.f
index eda8b61..bd3cfba 100644
--- a/src/dqagmv.f
+++ b/src/dqagmv.f
@@ -1,6 +1,7 @@
 c
 c
 c     Integration routine dqags.f from quadpack and dependencies: BEGIN
+c     Modified version for functions f(x,yi) with more than one variable
 c
 c
       subroutine dqagsmv(f,a,b,apar,np,epsabs,epsrel,result,abserr,
diff --git a/src/gray.f b/src/gray.f
index 68acc9c..5c9d3bc 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -881,7 +881,7 @@ c
 
       if (iprof.eq.1) then
         nprof=98
-        lprf=length(filenmprf)
+        lprf=len_trim(filenmprf)
         open(file=filenmprf(1:lprf)//'.prf',
      .                    status= 'unknown',unit=nprof)
         call profdata
@@ -892,7 +892,7 @@ c     read equilibrium data  from file if iequil=2
 c
       if (iequil.eq.2) then
         neqdsk=99
-        leqq=length(filenmeqq)
+        leqq=len_trim(filenmeqq)
          open(file=filenmeqq(1:leqq)//'.eqdsk',
      .        status= 'unknown', unit=neqdsk)
         call equidata
@@ -1081,7 +1081,7 @@ c
 c     initial beam data are measured in mm -> transformed to cm
 c
       nfbeam=97
-      lbm=length(filenmbm)
+      lbm=len_trim(filenmbm)
       open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam)
 c
       read(nfbeam,*) nisteer
diff --git a/src/green_func_p.f90 b/src/green_func_p.f90
index e222928..60c6395 100644
--- a/src/green_func_p.f90
+++ b/src/green_func_p.f90
@@ -8,7 +8,7 @@
 ! the current drive value by adjoint approach
 !
 !########################################################################
- USE const_and_precisions
+ USE const_and_precisions, only : wp_
 !-------
  IMPLICIT NONE
  CHARACTER(Len=1), PRIVATE :: adj_appr(6)    ! adjoint approach switcher
@@ -170,6 +170,7 @@
 !   K   - Spitzer's function
 !  dKdu = dK/du, i.e. its derivative over normalized momentum
 !=======================================================================
+ use const_and_precisions, only : comp_eps
  IMPLICIT NONE
  REAL(wp_), INTENT(in)  :: mu,Zeff,fc,u,q,gam
  REAL(wp_), INTENT(out) :: K,dKdu
@@ -234,6 +235,7 @@
 ! sfd(1),...,sfd(4) - coefficients of the polynomial expansion of the 
 !                 "Spitzer"-function (the same as in the Hirshman paper)
 !=======================================================================
+ use const_and_precisions, only : mc2_
  IMPLICIT NONE
  REAL(wp_), INTENT(in) :: mu,Zeff,fc
  INTEGER   :: n,i,j
@@ -385,6 +387,7 @@
 !  K   - Spitzer function
 ! dKdu - its derivative
 !=======================================================================
+ use const_and_precisions, only : zero,one
  IMPLICIT NONE
  REAL(wp_), INTENT(in)  :: Zeff,fc,u,q,gam
  REAL(wp_), INTENT(out) :: K,dKdu