diff --git a/Makefile.jetto b/Makefile.jetto
new file mode 100644
index 0000000..2866f26
--- /dev/null
+++ b/Makefile.jetto
@@ -0,0 +1,136 @@
+# Makefile for JETTO fortran 90 library 'gray'
+# Author: Lorenzo Figini (figini@ifp.cnr.it)
+#
+# Derived from
+# Makefile for JETTO fortran 77 library 'frantic'
+# Author: Derek Harting (d.harting@fz-juelich.de)
+# and
+# Makefile for ITCequ Fortran 90 library - definitions
+# P. Strand, elfps@chalmers.se
+# G. Corrigan, gcor@jet.uk
+# D. Harting, d.harting@fz-juelich.de
+# Set the environment from the top-level Makefile of JETTO
+# --------------------------------------------------------
+#include ../include.mk
+
+# Alternative search paths
+vpath %.f90 src
+
+JLIBDIR = ./
+F90 = gfortran
+DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
+
+# library name
+# ------------
+LIBNAME=$(JLIBDIR)/libgray.a
+
+# List source and object files
+# ----------------------------
+FSRC=$(wildcard src/*.f90)
+FOBJ0=$(subst src/,,$(FSRC))
+FOBJ=$(FOBJ0:.f90=.o)
+
+########################################
+# Set up compiler specific environment #
+########################################
+# set default compiler
+# --------------------
+FC=$(F90)
+
+# Set compiler debug and optimization flags
+# -----------------------------------------
+ifeq ("$(DBG)","")
+#     FFLAGS= $(AUTO) -O3 -Mpreprocess
+     FFLAGS=-O3 -cpp
+else
+#     FFLAGS= $(AUTO) -O0 -Mpreprocess -g -Minform=inform -Mbounds -Mchkptr
+     FFLAGS=-Wall -g -finit-real=nan -ffpe-trap=invalid -fcheck=all -fbounds-check
+endif
+
+# Set include directories
+# -----------------------
+#INC=-I$(JCOMMON_STD)
+INC=
+
+#####################################
+# Set up RULES specific environment #
+#####################################
+
+# suffixes
+# --------
+.SUFFIXES: .f90 .o .mod
+
+# PHONY targets
+# -------------
+.PHONY: all clean realclean $(MASTER_RULES) gray
+
+# default target
+# --------------
+all:
+	@($(MAKE) -C ../ all) || exit $$?
+
+# Set rules passed to top-level Makefile
+# --------------------------------------
+$(MASTER_RULES):
+	@($(MAKE) -C ../ $@) || exit $$?
+
+# rule for libgray.a
+# --------------------
+$(LIBNAME): $(FOBJ)
+	ar vr $(LIBNAME) $?
+
+# synonym for libgray.a
+# -----------------------
+gray: $(LIBNAME)
+
+# create rule for compiling f90 files
+# -----------------------------------
+%.o: %.f90
+	$(FC) -c $(DIRECTIVES) $(FFLAGS) $(INC)  $<
+
+# make 'realclean' option
+# -----------------------
+realclean:
+	rm -f *.o *.mod $(LIBNAME)
+
+# make 'clean' option
+# -------------------
+clean:
+	rm -f *.o
+
+# Dependencies
+# ------------
+main.o: const_and_precisions.o beams.o coreprofiles.o equilibrium.o \
+        graycore.o gray_params.o reflections.o
+gray_jetto1beam.o: const_and_precisions.o beams.o graycore.o gray_params.o \
+                   units.o
+
+graycore.o: const_and_precisions.o beamdata.o beams.o coreprofiles.o \
+            dispersion.o eccd.o equilibrium.o errcodes.o gray_params.o \
+            pec.o polarization.o limiter.o units.o utils.o
+beams.o: const_and_precisions.o dierckx.o reflections.o simplespline.o utils.o
+beamdata.o: const_and_precisions.o gray_params.o 
+conical.o: const_and_precisions.o
+coreprofiles.o: const_and_precisions.o dierckx.o gray_params.o simplespline.o \
+                utils.o
+dierckx.o: const_and_precisions.o
+dispersion.o: const_and_precisions.o eierf.o errcodes.o math.o quadpack.o
+eccd.o: const_and_precisions.o conical.o dierckx.o errcodes.o magsurf_data.o \
+        numint.o
+eierf.o: const_and_precisions.o
+errcodes.o: const_and_precisions.o
+gray_params.o: const_and_precisions.o utils.o
+equilibrium.o: const_and_precisions.o dierckx.o limiter.o minpack.o \
+               reflections.o simplespline.o utils.o gray_params.o
+magsurf_data.o: const_and_precisions.o gray_params.o equilibrium.o dierckx.o \
+                reflections.o simplespline.o units.o utils.o
+math.o: const_and_precisions.o
+minpack.o: const_and_precisions.o
+numint.o: const_and_precisions.o
+pec.o: const_and_precisions.o beamdata.o equilibrium.o gray_params.o \
+       magsurf_data.o utils.o
+polarization.o: const_and_precisions.o
+quadpack.o: const_and_precisions.o
+reflections.o: const_and_precisions.o limiter.o utils.o
+simplespline.o: const_and_precisions.o
+utils.o: const_and_precisions.o
diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90
new file mode 100644
index 0000000..c0f0d4e
--- /dev/null
+++ b/src/gray_jetto1beam.f90
@@ -0,0 +1,86 @@
+subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
+    nbnd, rbnd, zbnd, nrho, psrad, fpol, te, dne, zeff, qpsi, ibeam,   &
+    powin, alphain, betain, dpdv, jcd, pabs, icd, ierr)
+  use const_and_precisions, only : wp_
+  use units, only : ucenr,usumm,uprj0,uwbm,udisp,ubres,ucnt,uoutr,ueq,uprfin,  &
+      uflx,upec,uprm,ubeam
+  use gray_params, only : read_params,rtrparam_type,hcdparam_type,antctrl_type,&
+      eqparam_type,prfparam_type,outparam_type
+  use beams, only : read_beam2
+  use graycore, only : gray_main
+  implicit none
+! arguments
+  integer, intent(in) :: ijetto, mr, mz, nbnd, nrho, ibeam
+  real(wp_), dimension(mr), intent(in) :: r
+  real(wp_), dimension(mz), intent(in) :: z
+  real(wp_), dimension(mr,mz), intent(in) :: psin
+  real(wp_), intent(in) :: psia, rax, zax, powin, alphain, betain
+  real(wp_), dimension(nbnd), intent(in) :: rbnd, zbnd
+  real(wp_), dimension(nrho), intent(in) :: psrad, fpol, te, dne, zeff, qpsi
+  real(wp_), dimension(nrho), intent(out) :: dpdv, jcd
+  real(wp_), intent(out) :: pabs, icd
+  integer, intent(out) :: ierr
+! local variables
+  real(wp_), dimension(nrho) :: psinr
+  integer :: iox0
+  real(wp_) :: r0m,rvac,alpha0,beta0,psipol0,chipol0,p0mw,rwallm
+  real(wp_) :: fghz,x0,y0,z0,w1,w2,ri1,ri2,phiw,phir
+  real(wp_), dimension(5) :: rlim,zlim
+  logical, save :: firstcall=.true.
+  type(rtrparam_type) :: rtrp
+  type(hcdparam_type) :: hcdp
+  type(antctrl_type) :: antp
+  type(eqparam_type) :: eqp
+  type(prfparam_type) :: prfp
+  type(outparam_type) :: outp
+
+! if first call read parameters from external file
+  if (firstcall) then
+    call read_params('gray.data',rtrp,hcdp,antp,eqp,rwallm,prfp,outp,uprm)
+    rwallm=r(1)
+    antp%filenm='graybeam.data'
+    eqp%filenm='JETTO'
+    eqp%iequil=ijetto+1
+    prfp%filenm='JETTO'
+    outp%ipec=1
+    firstcall=.false.
+  end if
+  
+  rvac=rax
+  psinr=psrad
+
+  alpha0=alphain
+  beta0=betain
+  call read_beam2(antp%filenm,ibeam,alpha0,beta0,fghz,antp%iox,x0,y0,z0, &
+      w1,w2,ri1,ri2,phiw,phir,ubeam)
+  psipol0=antp%psi
+  chipol0=antp%chi
+  iox0=antp%iox
+  p0mw=powin*1.e-6_wp_
+
+! set simple limiter
+  r0m=sqrt(x0**2+y0**2)*0.01_wp_
+  call range2rect(rwallm,max(r0m,r(mr)),z(1),z(mz),rlim,zlim)
+
+! call main subroutine for the ibeam-th beam
+  call gray_main(r,z,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd,eqp,     &
+      psrad,te,dne,zeff,prfp, rlim,zlim, p0mw,fghz,alpha0,beta0,(/x0,y0,z0/),  &
+      w1,w2,ri1,ri2,phiw,phir,iox0,psipol0,chipol0, dpdv,jcd,pabs,icd,outp,    &
+      rtrp,hcdp,ierr,sqrt(psrad))
+
+! close output (debug) files
+  close(ucenr)
+  close(usumm)
+  close(uprj0)
+  close(uprj0+1)
+  close(uwbm)
+  close(udisp)
+  close(ubres)
+  close(ucnt)
+  close(uoutr)
+  close(ueq)
+  close(uprfin)
+  close(uflx)
+  close(upec)
+
+end subroutine gray_jetto1beam
\ No newline at end of file
diff --git a/src/graycore.f90 b/src/graycore.f90
index d7ee8f5..7349ff3 100644
--- a/src/graycore.f90
+++ b/src/graycore.f90
@@ -7,7 +7,7 @@ contains
   subroutine gray_main(rv,zv,psin,psia,psinr,fpol,qpsi,rvac,rax,zax,rbnd,zbnd, &
               eqp,psrad,terad,derad,zfc,prfp,  rlim,zlim,                      &
               p0,fghz,alpha0,beta0,xv0,w1,w2,ri1,ri2,phiw,phir,iox0,           &
-              psipol0,chipol0,  dpdv,jcd,pabs,icd,  outp,rtrp,hcdp,ierr)
+              psipol0,chipol0,  dpdv,jcd,pabs,icd,  outp,rtrp,hcdp,ierr, rhout)
     use const_and_precisions, only : zero, one
     use coreprofiles, only : set_prfan, set_prfspl, temp, fzeff
     use dispersion, only : expinit
@@ -33,10 +33,10 @@ contains
     type(rtrparam_type), intent(in) :: rtrp
     type(hcdparam_type), intent(in) :: hcdp
 
-    real(wp_), dimension(:), allocatable, intent(in) :: psrad, terad, derad, zfc
-    real(wp_), dimension(:), allocatable, intent(in) :: rv, zv, psinr, fpol, qpsi
-    real(wp_), dimension(:), allocatable, intent(in) :: rbnd, zbnd, rlim, zlim
-    real(wp_), dimension(:,:), allocatable, intent(in) :: psin
+    real(wp_), dimension(:), intent(in) :: psrad, terad, derad, zfc
+    real(wp_), dimension(:), intent(in) :: rv, zv, psinr, fpol, qpsi
+    real(wp_), dimension(:), intent(in) :: rbnd, zbnd, rlim, zlim
+    real(wp_), dimension(:,:), intent(in) :: psin
     real(wp_), intent(in) :: psia, rvac, rax, zax
     integer, intent(in) :: iox0
     real(wp_), intent(in) :: p0, fghz, psipol0, chipol0
@@ -45,6 +45,7 @@ contains
 
     real(wp_), intent(out) :: pabs,icd
     real(wp_), dimension(:), intent(out) :: dpdv,jcd
+    real(wp_), dimension(:), intent(in), optional :: rhout
 
     integer, intent(out) :: ierr
 ! local variables
@@ -236,7 +237,7 @@ contains
     if(nray > 1) call print_projxyzt(st,yw,1)
 
 ! compute power and current density profiles for all rays
-    call pec_init(ipec) !,sqrt(psinr))
+    call pec_init(ipec,rhout)
     nnd=size(rhop_tab)
     allocate(jphi(nnd),pins(nnd),currins(nnd))
     call spec(psjki,ppabs,ccci,iiv,pabs,icd,dpdv,jphi,jcd,pins,currins)
@@ -264,7 +265,7 @@ contains
 !    call dealloc_beam(yw,ypw,xc,du1,gri,ggri,psjki,ppabs,ccci, &
 !                      tau0,alphaabs0,dids0,ccci0,p0jk,ext,eyt,iiv)
 !    call dealloc_pec
-!    deallocate(jphi,pins,currins)
+    deallocate(jphi,pins,currins)
 ! ======= free memory  END  ======
   end subroutine gray_main
 
diff --git a/src/units.f90 b/src/units.f90
index d2db0a0..90580bd 100644
--- a/src/units.f90
+++ b/src/units.f90
@@ -1,13 +1,13 @@
 module units
 ! STANDARD
-  integer, parameter :: uprm   =  2, ubeam = 97, uprf  = 98, ueq    = 99
-  integer, parameter :: uprfin = 55, uflx  = 56, ubres = 70, ucnt   = 71
-  integer, parameter :: uprj0  =  8, uwbm  = 12, ucenr =  4, uoutr  = 33
-  integer, parameter :: udisp  = 17, upec  = 48, usumm =  7
+!  integer, parameter :: uprm   =  2, ubeam = 97, uprf  = 98, ueq    = 99
+!  integer, parameter :: uprfin = 55, uflx  = 56, ubres = 70, ucnt   = 71
+!  integer, parameter :: uprj0  =  8, uwbm  = 12, ucenr =  4, uoutr  = 33
+!  integer, parameter :: udisp  = 17, upec  = 48, usumm =  7
 ! JINTRAC
-!  integer, parameter :: uprm   =602, ubeam =603, uprf  =644, ueq    =644
-!  integer, parameter :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
-!  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
-!  integer, parameter :: udisp  =617, upec  =648, usumm =607, uhead  =638
+  integer, parameter :: uprm   =602, ubeam =603, uprf  =644, ueq    =644
+  integer, parameter :: uprfin =645, uflx  =646, ubres =630, ucnt   =631
+  integer, parameter :: uprj0  =608, uwbm  =612, ucenr =604, uoutr  =633
+  integer, parameter :: udisp  =617, upec  =648, usumm =607, uhead  =638
 
 end module units
\ No newline at end of file
diff --git a/srcjetto/gray.f b/srcjetto/gray.f
index 84cce7a..a2f1611 100644
--- a/srcjetto/gray.f
+++ b/srcjetto/gray.f
@@ -2,11 +2,6 @@
       subroutine gray(ijetto, mr, mz, mrd, r, z, psin, psiax, psibnd,
      .    rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff,
      .    qsf, nbeam, powin, alphin, betain, dpdv, jcd, pec, icd, ierr)
-      use units, only : ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres,
-     .    ucnt,uoutr,ueq,uprfin,uflx,upec
-      use gray_params, only : read_params
-      use beams, only : read_beam2
-      use graycore, only : gray_main
       implicit none
 !   input arguments
       integer ijetto, mr, mz, nbnd, nrho, nbeam
@@ -21,10 +16,6 @@
 !   gray_main output arguments
       real*8 dpdvloop(nrho), jcdloop(nrho), pecloop, icdloop
       integer ierr
-!   local variables
-      real*8 rlim(5),zlim(5)
-      logical firstcall=.true.
-      save firstcall
 
 ! === input arguments ==================================================
 !
@@ -86,17 +77,6 @@
 !
 ! ======================================================================
 
-!   if first call read parameters from external file
-      if (firstcall) then
-        call read_params('gray.data',rtrp,hcdp,antp,eqp,rwallm,
-     .    prfp,outp,uprm)
-        antp%filenm='graybeam.data'
-        eqp%filenm='JETTO'
-        eqp%iequil=ijetto+1
-        prfp%filenm='JETTO'
-        firstcall=.false.
-      end if
-
 !   set output variables to 0
       do i=1,nrho
         dpdv(i) = 0.d0
@@ -111,26 +91,12 @@
 
 !       read j-th beam properties from file
 !       and adjust alpha/beta if out of the allowed range
-        alpha0=alphin(j)
-        beta0=betain(j)
-        p0mw=powin(j)*1.d-6
-        call read_beam2(antp%filenm,j,alpha0,beta0,fghz,antp%iox,
-     .                  x0,y0,z0,w1,w2,ri1,ri2,phiw,phir,ubeam)
-        psipol0=antp%psi
-        chipol0=antp%chi
-        iox0=antp%iox
-
-!       set simple limiter
-        r0m=sqrt(x0**2+y0**2)*0.01d0
-        call range2rect(rwallm,max(r0m,rv(mr)),zv(1),zv(mz),rlim,zlim)
 
 !       call main subroutine for the j-th beam
-        subroutine gray_main(r,z,psin(1:mr,:),psibnd-psiax,
-     .      psijet,-f,-qsf,rax,rax,zax,rbnd,zbnd,eqp,
-     .      psijet,te,dne,zeff,prfp,rlim,zlim,
-     .      p0mw,fghz,alpha0,beta0,(/x0,y0,z0/),
-     .      w1,w2,ri1,ri2,phiw,phir,iox0,psipol0,chipol0,
-     .      dpdvloop,jcdloop,pecpool,icdloop,outp,rtrp,hcdp,ierr)
+        subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin(1:mr,:),
+     .    psibnd-psiax, rax, zax, nbnd, rbnd, zbnd, nrho, psijet, -f,
+     .    te, dne, zeff, -qsf, j, powin(j), alphin(j), betain(j),
+     .    dpdvloop, jcdloop, pecloop, icdloop, ierr)
 
 !       add contribution of j-th beam to the total
 !       adapting output data to JETTO convention on toroidal angle
@@ -144,9 +110,5 @@
 !     end of loop over beams with power>0
       end do
 
-!     close output (debug) files
-      close(ucenr,usumm,uprj0,uprj0+1,uwbm,udisp,ubres,ucnt,uoutr,
-     .    ueq,uprfin,uflx,upec)
-
       return
       end