diff --git a/Makefile b/Makefile
index 97ebf1a..077048a 100644
--- a/Makefile
+++ b/Makefile
@@ -6,7 +6,7 @@ MAINOBJ=gray.o
 OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \
     graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
     interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
-    reflections.o simplespline.o utils.o
+    reflections.o simplespline.o utils.o beamdata.o
     
 # Alternative search paths
 vpath %.f90 src
@@ -28,7 +28,7 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
 gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \
         graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
         interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
-        reflections.o simplespline.o utils.o
+        reflections.o simplespline.o utils.o beamdata.o
 conical.o: const_and_precisions.o
 dierckx.o: const_and_precisions.o
 dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
@@ -45,6 +45,7 @@ quadpack.o: const_and_precisions.o
 reflections.o: const_and_precisions.o utils.o
 simplespline.o: const_and_precisions.o
 utils.o: const_and_precisions.o
+beamdata.o: const_and_precisions.o
 
 # General object compilation command
 %.o: %.f90
diff --git a/src/beamdata.f90 b/src/beamdata.f90
new file mode 100644
index 0000000..03831b5
--- /dev/null
+++ b/src/beamdata.f90
@@ -0,0 +1,80 @@
+module beamdata
+  use const_and_precisions, only : wp_
+  implicit none
+
+  integer, parameter :: jmx=31,kmx=36,nmx=8000
+  integer, save :: nrayr,nrayth,nstep
+  real(wp_), dimension(:,:,:), allocatable, save :: psjki,ppabs,ccci,tauv,alphav
+  real(wp_), dimension(:,:,:), allocatable, save :: pdjki,currj,didst
+  integer, dimension(:,:), allocatable, save :: iiv,iop,iow,ihcd,istore
+  real(wp_), dimension(:,:), allocatable, save :: tau1v
+  real(wp_), dimension(:), allocatable, save :: q
+  real(wp_), dimension(:,:,:), allocatable, save :: yyrfl !(:,:,6)
+  real(wp_), dimension(:,:,:), allocatable, save :: ywrk,ypwrk !(6,:,:)
+  real(wp_), dimension(:,:,:), allocatable, save :: xc,xco,du1,du1o !(3,:,:)
+  real(wp_), dimension(:,:,:), allocatable, save :: gri,dgrad2v !(3,:,:)
+  real(wp_), dimension(:,:,:,:), allocatable, save :: ggri !(3,3,:,:)
+  real(wp_), dimension(:,:), allocatable, save :: grad2
+  real(wp_), dimension(:), allocatable, save :: dffiu,ddffiu
+  complex(wp_), dimension(:,:,:), allocatable, save :: ext,eyt
+
+contains
+
+  subroutine alloc_beam(ierr)
+    implicit none
+    integer, intent(out) :: ierr
+    
+    call dealloc_beam
+    allocate(psjki(nrayr,nrayth,nstep),  ppabs(nrayr,nrayth,nstep), &
+             pdjki(nrayr,nrayth,nstep),   ccci(nrayr,nrayth,nstep), &
+             currj(nrayr,nrayth,nstep),  didst(nrayr,nrayth,nstep), &
+              tauv(nrayr,nrayth,nstep), alphav(nrayr,nrayth,nstep), &
+                     iiv(nrayr,nrayth),          iop(nrayr,nrayth), &
+                     iow(nrayr,nrayth),        tau1v(nrayr,nrayth), &
+                    ihcd(nrayr,nrayth),       istore(nrayr,nrayth), &
+                              q(nrayr),      yyrfl(nrayr,nrayth,6), &
+                  ywrk(6,nrayr,nrayth),      ypwrk(6,nrayr,nrayth), &
+                    xc(3,nrayr,nrayth),        xco(3,nrayr,nrayth), &
+                   du1(3,nrayr,nrayth),       du1o(3,nrayr,nrayth), &
+                   gri(3,nrayr,nrayth),    dgrad2v(3,nrayr,nrayth), &
+                ggri(3,3,nrayr,nrayth),        grad2(nrayr,nrayth), &
+                          dffiu(nrayr),              ddffiu(nrayr), &
+                 ext(nrayr,nrayth,0:4),      eyt(nrayr,nrayth,0:4), &
+                             stat=ierr)
+    if (ierr/=0) call dealloc_beam
+  end subroutine alloc_beam
+
+  subroutine dealloc_beam
+    implicit none
+    if (allocated(psjki))   deallocate(psjki)
+    if (allocated(ppabs))   deallocate(ppabs)
+    if (allocated(pdjki))   deallocate(pdjki)
+    if (allocated(ccci))    deallocate(ccci)
+    if (allocated(currj))   deallocate(currj)
+    if (allocated(didst))   deallocate(didst)
+    if (allocated(tauv))    deallocate(tauv)
+    if (allocated(alphav))  deallocate(alphav)
+    if (allocated(iiv))     deallocate(iiv)
+    if (allocated(iop))     deallocate(iop)
+    if (allocated(iow))     deallocate(iow)
+    if (allocated(ihcd))    deallocate(ihcd)
+    if (allocated(istore))  deallocate(istore)
+    if (allocated(tau1v))   deallocate(tau1v)
+    if (allocated(q))       deallocate(q)
+    if (allocated(yyrfl))   deallocate(yyrfl)
+    if (allocated(ywrk))    deallocate(ywrk)
+    if (allocated(ypwrk))   deallocate(ypwrk)
+    if (allocated(xc))      deallocate(xc)
+    if (allocated(xco))     deallocate(xco)
+    if (allocated(du1))     deallocate(du1)
+    if (allocated(du1o))    deallocate(du1o)
+    if (allocated(gri))     deallocate(gri)
+    if (allocated(dgrad2v)) deallocate(dgrad2v)
+    if (allocated(ggri))    deallocate(ggri)
+    if (allocated(grad2))   deallocate(grad2)
+    if (allocated(dffiu))   deallocate(dffiu)
+    if (allocated(ddffiu))  deallocate(ddffiu)
+    if (allocated(ext))     deallocate(ext)
+    if (allocated(eyt))     deallocate(eyt)
+  end subroutine dealloc_beam
+end module beamdata
diff --git a/src/gray.f b/src/gray.f
index ceb2d57..e490d01 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -66,6 +66,7 @@ c
         pabstott=pabstott+pabstot
         currtott=currtott+currtot
       end if
+      call vectfree
       print*,' '
       write(6,*) 'Pabs (MW), Icd (kA) = ', pabstott,currtott*1.0e3_wp_
 c
@@ -76,17 +77,17 @@ c
       subroutine gray_integration
       use const_and_precisions, only : wp_
       use graydata_flags, only : dst
+      use beamdata, only : nstep
       implicit none
 c local variables
       integer :: i
       real(wp_) :: st0
 c common/external functions/variables
-      integer :: istep,nstep,istop,index_rt
+      integer :: istep,istop,index_rt
       real(wp_) :: st,strfl11
 c
       common/ss/st
       common/istep/istep
-      common/nstep/nstep
       common/istop/istop
       common/strfl11/strfl11
       common/index_rt/index_rt
@@ -119,16 +120,16 @@ c
       use graydata_flags, only : ibeam,iwarm,iequil,iprof,
      . filenmeqq,filenmprf,filenmbm
       use graydata_anequil, only : dens0,te0
+      use beamdata, only : nrayr
       implicit none
 c local variables
       integer :: iproj,nfilp
       real(wp_) :: currtka,pabs,currt
 c common/external functions/variables
-      integer :: nrayr,nrayth,index_rt
+      integer :: index_rt
       real(wp_) :: st,taumn,taumx,pabstot,currtot
 c
       common/ss/st
-      common/nray/nrayr,nrayth
       common/index_rt/index_rt
       common/taumnx/taumn,taumx,pabstot,currtot
 c
@@ -180,46 +181,34 @@ c
       use graydata_flags, only : iwarm,istpr0,istpl0,dst,ipass,igrad
       use graydata_par, only : psipol0,chipol0
       use interp_eqprof, only : zbmin,zbmax,rlim,zlim,nlim
+      use beamdata, only : nrayr,nrayth,psjki,ppabs,ccci,iiv,tauv,ihcd,
+     . istore,iop,iow,tau1v,yyrfl,ywrk,ypwrk,ext,eyt
       implicit none
 c local constants
-      integer, parameter :: jmx=31,kmx=36,nmx=8000
       real(wp_), parameter :: taucr=12.0_wp_
 c arguments
       integer :: i, istop
 c local variables
-      integer iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1,k1,
-     . kkkk
+      integer :: iproj,nfilp,irfl,j,k,kkk,iopmin,iowmax,iowmin,ivac,j1,
+     . k1,kkkk
       real(wp_) :: qqout,uuout,vvout,qqin2,uuin2,vvin2,rrm,zzm,anvjkr,
      . aknmin,akdotn,anwclr
       real(wp_), dimension(6) :: y,dery
       real(wp_), dimension(3) :: xvrfl,anvrfl,xvvac,anw
-      real(wp_), dimension(3,jmx,kmx) :: xvjk,anvjk
+      real(wp_), dimension(3,nrayr,nrayth) :: xvjk,anvjk
       complex(wp_) :: extr,eytr,exin2,eyin2
       logical :: ins_pl
 c common/external functions/variables
-      integer :: nrayr,nrayth,istpr,istpl,jclosest,ierr,index_rt
-      integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore
+      integer :: istpr,istpl,jclosest,ierr,index_rt
       real(wp_) :: psinv,psinv11,taumn,taumx,pabstot,currtot,psipol,
      . chipol,powrfl,strfl11
       real(wp_), dimension(3) :: anwcl,xwcl,xv,anv
-      real(wp_), dimension(jmx,kmx) :: tau1v
-      real(wp_), dimension(jmx,kmx,6) :: yyrfl
-      real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk
-      real(wp_), dimension(jmx,kmx,nmx) :: ppabs,ccci,tauv,alphav,psjki
-      complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt
       logical :: inside_plasma
 c
       external inside_plasma
 c
-      common/pcjki/ppabs,ccci
-      common/atjki/tauv,alphav
-      common/tau1v/tau1v
-      common/nray/nrayr,nrayth
       common/istgr/istpr,istpl
-      common/iiv/iiv
-      common/iov/iop,iow,ihcd,istore
       common/refln/anwcl,xwcl,jclosest
-      common/psjki/psjki
       common/psival/psinv
       common/psinv11/psinv11
       common/ierr/ierr
@@ -227,12 +216,9 @@ c
       common/xv/xv
       common/anv/anv
       common/polcof/psipol,chipol
-      common/evt/ext,eyt
-      common/yyrfl/yyrfl
       common/powrfl/powrfl
       common/strfl11/strfl11
       common/index_rt/index_rt
-      common/wrk/ywrk,ypwrk
 c
       pabstot=0.0_wp_
       currtot=0.0_wp_
@@ -476,41 +462,30 @@ c
       use const_and_precisions, only : wp_,pi
       use graydata_flags, only : istpl0
       use graydata_par, only : psdbnd
+      use beamdata, only : ywrk,psjki,tauv,alphav,pdjki,
+     . currj,didst,q,tau1v,nrayr,nrayth
       implicit none
 c local constants
-      integer, parameter :: ndim=6,jmx=31,kmx=36,nmx=8000
+      integer, parameter :: ndim=6
       real(wp_), parameter :: taucr=12.0_wp_
 c arguments
-      integer i,j,k
+      integer :: i,j,k
 c local variables
       real(wp_) :: x,y,z,rr,rrm,zzm,stm,xxm,yym,anx,any,anz,anr,anphi,
      . phi,phideg,psi,rhot,bbr,bbz,bpol,dens11,dids11,dpdv11,ajcd11,
      . alpha11,taujk,tau11,pt11
 c common/external functions/variables
-      integer :: nharm,nhf,iohkawa,index_rt,nrayr,nrayth,istpr,istpl
+      integer :: nharm,nhf,iohkawa,index_rt,istpr,istpl
       real(wp_) :: p0mw,st,brr,bphi,bzz,ajphi,btot,xg,yg,dens,ddens,
      . tekev,alpha,effjcd,akim,tau0,anpl,anpr,ddr,an2s,an2,fdia,bdotgr,
      . ddi,ddr11,anprr,anpri,frhotor
-      real(wp_), dimension(jmx) :: q
-      real(wp_), dimension(jmx,kmx) :: tau1v
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk
-      real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki,
-     . ppabs,currj,didst,ccci
       complex(wp_) :: ex,ey,ez
 c
-      common/psjki/psjki
-      common/atjki/tauv,alphav
-      common/dpjjki/pdjki,currj
-      common/pcjki/ppabs,ccci
-      common/dcjki/didst
       common/nharm/nharm,nhf
       common/iokh/iohkawa
       common/p0/p0mw
-      common/tau1v/tau1v
-      common/qw/q
       common/index_rt/index_rt
       common/ss/st
-      common/nray/nrayr,nrayth
       common/istgr/istpr,istpl
       common/parpl/brr,bphi,bzz,ajphi
       common/btot/btot
@@ -523,7 +498,6 @@ c
       common/nplr/anpl,anpr
       common/ddd/ddr,an2s,an2,fdia,bdotgr,ddi,ddr11
       common/nprw/anprr,anpri
-      common/wrk/ywrk,ypwrk
 c
       x=ywrk(1,j,k)
       y=ywrk(2,j,k)
@@ -660,9 +634,8 @@ c
       use graydata_anequil
       use interp_eqprof, only : rmxm,rlim,zlim,nlim,zbmin,zbmax,
      . btrcen,rcen,alloc_lim
+      use beamdata, only : nrayr,nrayth,nstep
       implicit none
-c local constants
-      integer, parameter :: nmx=8000
 c local variables
       integer :: nfil,iox,ierr,leqq,lprf
       real(wp_) :: dummy,bresg,r00,anr0c,anphi0c,fghz,zeff,
@@ -670,13 +643,10 @@ c local variables
       character(len=8) :: wdat
       character(len=10) :: wtim
 c common/external functions/variables
-      integer :: nstep,nrayr,nrayth
       real(wp_) :: xgcn,ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw,
      . phir,anx0c,any0c,anz0c,x00,y00,z00,bres,p0mw,sox,alpha0,beta0
 c
       common/xgcn/xgcn
-      common/nstep/nstep
-      common/nray/nrayr,nrayth
       common/parwv/ak0,akinv,fhz
       common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
       common/anic/anx0c,any0c,anz0c
@@ -946,9 +916,9 @@ c     set simple limiter as two cylindrical walls at rwallm and r00
         rlim(3)=rlim(2)
         rlim(4)=rlim(1)
         rlim(5)=rlim(1)
-        zlim(1)=(z00-dst*nmx)*1.0e-2_wp_
+        zlim(1)=(z00-dst*nstep)*1.0e-2_wp_
         zlim(2)=zlim(1)
-        zlim(3)=(z00+dst*nmx)*1.0e-2_wp_
+        zlim(3)=(z00+dst*nstep)*1.0e-2_wp_
         zlim(4)=zlim(3)
         zlim(5)=zlim(1)
         ipass=abs(ipass)
@@ -1906,9 +1876,9 @@ c
       use utils, only : locate, intlin
       implicit none
 c arguments
-      real(wp_) qval,psival
+      real(wp_) :: qval,psival
 c local variables
-      integer ncnt,i1,ipr
+      integer :: ncnt,i1,ipr
       real(wp_), dimension(npoints) :: rcn,zcn
 c
       ncnt=(npoints-1)/2
@@ -1933,9 +1903,9 @@ c
 c local constants
       integer, parameter :: icmx=2002
 c local variables
-      integer j,k,n,nconts,inc,nctot
+      integer :: j,k,n,nconts,inc,nctot
       integer, dimension(10) :: ncpts
-      real(wp_) btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
+      real(wp_) :: btmx,btmn,zzk,rrj,bbphi,bbr,bbz,bbb
       real(wp_), dimension(icmx) :: rrcb,zzcb
 c common/external functions/variables
       real(wp_) :: bres
@@ -1988,9 +1958,9 @@ c
       use dierckx, only : curfit,splev,splder
       implicit none
 c local variables
-      integer ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2
+      integer :: ierr,i,iopt,ier,kspl,npest,lwrkf,nu,nn,nn1,nn2
       integer, dimension(:), allocatable :: iwrkf
-      real(wp_) aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi,
+      real(wp_) :: aat,aan,ffact,psrad0,terad0,derad0,zfc0,psradi,
      . teradi,deradi,zfci,xb,xe,sspl,dnpp,ddnpp,d2dnpp,dpsb,fp
       real(wp_), dimension(:), allocatable :: wf,wrkf,wrkfd,densi,
      . ddensi,d2densi
@@ -2101,13 +2071,13 @@ c
       use simplespline, only : difcsn
       implicit none
 c arguments
-      integer nnr
+      integer :: nnr
 c local variables
-      integer iopt,ier,k
-      real(wp_) dx,drhot
+      integer :: iopt,ier,k
+      real(wp_) :: dx,drhot
       real(wp_), dimension(nr) :: rhotnr
 c common/external functions/variables
-      real(wp_) rhotsx
+      real(wp_) :: rhotsx
 c
       common/rhotsx/rhotsx
 c
@@ -2148,10 +2118,10 @@ c
       use simplespline, only :spli
       implicit none
 c arguments
-      real(wp_) psi,fq_eq
+      real(wp_) :: psi,fq_eq
 c local variables
-      integer irt
-      real(wp_) dps
+      integer :: irt
+      real(wp_) :: dps
 c
       irt=int((nr-1)*psi+1)
       if(irt.eq.0) irt=1
@@ -2229,8 +2199,8 @@ c
 c
 c
       subroutine rhopol
-      use dierckx, only : curfit,splev
       use const_and_precisions, only : wp_
+      use dierckx, only : curfit,splev
       implicit none
 c local constants
       integer, parameter :: nnr=101,nrest=nnr+4,
@@ -2238,14 +2208,14 @@ c local constants
 c local variables
       integer :: i,iopt,ier,kspl
       integer, dimension(nrest) :: iwrkp
-      real(wp_) dr,psin,xb,xe,ss,rp
+      real(wp_) :: dr,psin,xb,xe,ss,rp
       real(wp_), dimension(nnr) :: rhop,rhot,rhopi
       real(wp_), dimension(lwrkp) :: wrkp
       real(wp_), dimension(nrest) :: wp
 c common/external functions/variables
       integer :: nsrp
-      real(wp_), dimension(nrest) :: trp,crp
       real(wp_) :: frhotor
+      real(wp_), dimension(nrest) :: trp,crp
 c
       common/coffrtp/trp
       common/coffrn/nsrp
@@ -2324,7 +2294,7 @@ c arguments
       real(wp_) :: h
       real(wp_), dimension(icmx) :: rcon,zcon
 c local variables
-      integer i,j,k,l,ico,nrqmax,nreq,nzeq,iclast,mpl,ix,jx,
+      integer :: i,j,k,l,ico,nrqmax,nreq,nzeq,iclast,mpl,ix,jx,
      . mxr,n1,jm,jfor,lda,ldb,jabs,jnb,kx,ikx,itm,inext,in
       integer, dimension(3,2) :: ja
       integer, dimension(1000) :: lx
@@ -2588,7 +2558,7 @@ c local constants
      . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54,
      . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
 c local variables
-      integer ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr
+      integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr
       integer, dimension(kwrk) :: iwrk
       real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator,
      . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height,height2,r2iav,currp,
@@ -2989,7 +2959,7 @@ c local constants
      . lwrk=4*(nnintp+nlam)+11*(njest+nlest)+njest*nnintp+nlest+54,
      . kwrk=nnintp+nlam+njest+nlest+3,lw01=nnintp*4+nlam*3+nnintp*nlam
 c local variables
-      integer ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr
+      integer :: ier,ierr,l,jp,ipr,jpr,inc,inc1,iopt,njp,nlm,ninpr
       integer, dimension(kwrk) :: iwrk
       real(wp_) :: lam,dlam,anorm,b2av,dvdpsi,dadpsi,ratio_cdator,
      . ratio_cdbtor,ratio_pltor,qq,fc,rhot2q,height2,r2iav,currp,area,
@@ -3359,8 +3329,8 @@ c
 c
 c
       subroutine rhopol_an
-      use dierckx, only : curfit,splev
       use const_and_precisions, only : wp_
+      use dierckx, only : curfit,splev
       use graydata_par, only : sgniphi
       use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam
       use interp_eqprof, only : psia
@@ -3438,8 +3408,8 @@ c
 c
 c
       function frhotor_an(rhop)
-      use dierckx, only : splev
       use const_and_precisions, only : wp_
+      use dierckx, only : splev
       implicit none
 c local contants
       integer, parameter :: nnr=101,nrest=nnr+4
@@ -3499,35 +3469,24 @@ c
       use const_and_precisions, only : wp_
       use dispersion, only: expinit
       use graydata_flags, only : iwarm
+      use beamdata, only : jmx,kmx,nmx,psjki,tauv,alphav,pdjki,ppabs,
+     . currj,didst,ccci,iiv,iop,iow,ihcd,istore,tau1v,nrayr,nrayth,
+     . nstep,alloc_beam
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36,nmx=8000
 c local variables
-      integer :: i,j,k
+      integer :: i,j,k,ierr
 c common/external functions/variables
-      integer nrayr,nrayth,nstep
-      integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore
       real(wp_) :: jclosest
       real(wp_), dimension(3) :: anwcl,xwcl
-      real(wp_), dimension(jmx,kmx) :: tau1v
-      real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki,
-     . ppabs,currj,didst,ccci
 c
-      common/iiv/iiv
-      common/iov/iop,iow,ihcd,istore
-      common/psjki/psjki
-      common/atjki/tauv,alphav
-      common/dpjjki/pdjki,currj
-      common/pcjki/ppabs,ccci
-      common/dcjki/didst
-      common/nray/nrayr,nrayth
-      common/nstep/nstep
-      common/tau1v/tau1v
       common/refln/anwcl,xwcl,jclosest
 c
       if(nstep.gt.nmx) nstep=nmx
       if(nrayr.gt.jmx) nrayr=jmx
       if(nrayth.gt.kmx) nrayth=kmx
+      call alloc_beam(ierr)
+      if (ierr.ne.0) return
+c
       jclosest=nrayr+1
       anwcl(1:3)=0.0_wp_
       xwcl(1:3)=0.0_wp_
@@ -3559,30 +3518,26 @@ c
       end
 c
 c
+c
+      subroutine vectfree
+      use beamdata, only : dealloc_beam
+      implicit none
+c
+      call dealloc_beam
+c
+      return
+      end
+c
+c
 c
       subroutine vectinit2
       use const_and_precisions, only : wp_
+      use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,
+     . didst,ccci,iiv,iop,iow,ihcd,nrayr,nrayth,nstep
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36,nmx=8000
 c local variables
       integer :: i,j,k
-c common/external functions/variables
-      integer nrayr,nrayth,nstep
-      integer, dimension(jmx,kmx) :: iiv,iop,iow,ihcd,istore
-      real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki,
-     . ppabs,currj,didst,ccci
 c
-      common/iiv/iiv
-      common/iov/iop,iow,ihcd,istore
-      common/psjki/psjki
-      common/atjki/tauv,alphav
-      common/dpjjki/pdjki,currj
-      common/pcjki/ppabs,ccci
-      common/dcjki/didst
-      common/nray/nrayr,nrayth
-      common/nstep/nstep
-c      
       do i=1,nstep
         do k=1,nrayth
           do j=1,nrayr
@@ -3611,7 +3566,7 @@ c
       use const_and_precisions, only : wp_
       implicit none
 c common/external functions/variables
-      integer istep,istpr,istpl,ierr,istop
+      integer :: istep,istpr,istpl,ierr,istop
 c
       common/istep/istep
       common/istgr/istpr,istpl
@@ -3631,20 +3586,10 @@ c
 c
       subroutine updatepos
       use const_and_precisions, only : wp_
+      use beamdata, only : xc,xco,du1,du1o,ywrk,nrayr,nrayth
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36
 c local variables
-      integer j,k
-c common/external functions/variables
-      integer :: nrayr,nrayth
-      real(wp_), dimension(3,jmx,kmx) :: xc,xco,du1,du1o
-      real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk
-c
-      common/nray/nrayr,nrayth
-      common/grco/xco,du1o
-      common/grc/xc,du1
-      common/wrk/ywrk,ypwrk
+      integer :: j,k
 c
       do j=1,nrayr
         do k=1,nrayth
@@ -3676,30 +3621,15 @@ c
 c
       subroutine gradi
       use const_and_precisions, only : wp_
+      use beamdata, only : dffiu,ddffiu,xc,xco,du1,du1o,gri,ggri,
+     . dgrad2v,grad2,nrayr,nrayth
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36
 c local variables
-      integer iv,j,jm,jp,k,km,kp
-      real(wp_) ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz,
+      integer :: iv,j,jm,jp,k,km,kp
+      real(wp_) :: ux,uxx,uxy,uxz,uy,uyy,uyz,uz,uzz,
      . dfu,dfuu,gx,gxx,gxy,gxz,gy,gyy,gyz,gz,gzz
       real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgu,
      . dgg1,dgg2,dgg3,df1,df2,df3
-c common/external functions/variables
-      integer :: nrayr,nrayth
-      real(wp_), dimension(jmx) :: dffiu,ddffiu
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: xc,xco,du1,du1o,gri,dgrad2v
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
-c
-      common/nray/nrayr,nrayth
-      common/fi/dffiu,ddffiu
-      common/grco/xco,du1o
-      common/grc/xc,du1
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
 c
 c     compute  grad u and grad(S_I)
 c
@@ -3851,7 +3781,7 @@ c
 c arguments
       real(wp_), dimension(3) :: dxv1,dxv2,dxv3,dgg
 c local variables
-      real(wp_) denom,aa1,aa2,aa3
+      real(wp_) :: denom,aa1,aa2,aa3
 c
       aa1=(dxv2(2)*dxv3(3)-dxv3(2)*dxv2(3))
       aa2=(dxv1(2)*dxv3(3)-dxv1(3)*dxv3(2))
@@ -3907,29 +3837,20 @@ c     Runge-Kutta integrator
 c
       use const_and_precisions, only : wp_
       use graydata_flags, only : dst,igrad
+      use beamdata, only : nrayr,nrayth,ywrk,ypwrk,grad2,dgrad2v,
+     . gri,ggri
       implicit none
 c parameter
-      integer, parameter :: ndim=6,jmx=31,kmx=36
+      integer, parameter :: ndim=6
 c local variables
-      integer ieq,iv,j,jv,k,kkk
-      real(wp_) h,hh,h6
+      integer :: ieq,iv,j,jv,k,kkk
+      real(wp_) :: h,hh,h6
       real(wp_), dimension(ndim) :: y,yy,fk1,fk2,fk3,fk4
 c common/external functions/variables
-      integer :: nrayr,nrayth
-      real(wp_) gr2
+      real(wp_) :: gr2
       real(wp_), dimension(3) :: dgr2,dgr
       real(wp_), dimension(3,3) :: ddgr
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: dgrad2v,gri
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
 c
-      common/nray/nrayr,nrayth
-      common/wrk/ywrk,ypwrk
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
       common/gr/gr2
       common/dgr/dgr2,dgr,ddgr
 c
@@ -3985,11 +3906,12 @@ c
       subroutine gwork(j,k)
       use const_and_precisions, only : wp_
       use graydata_flags, only : igrad
+      use beamdata, only : ywrk,ypwrk,grad2,dgrad2v,gri,ggri
       implicit none
 c local constants
-      integer, parameter :: ndim=6,jmx=31,kmx=36
+      integer, parameter :: ndim=6
 c arguments
-      integer j,k
+      integer :: j,k
 c local variables
       integer :: ieq,iv,jv
       real(wp_), dimension(ndim) :: yy,yyp
@@ -3997,16 +3919,7 @@ c common/external functions/variables
       real(wp_) :: gr2
       real(wp_), dimension(3) :: dgr2,dgr
       real(wp_), dimension(3,3) :: ddgr
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: dgrad2v,gri
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk,ypwrk
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
 c
-      common/wrk/ywrk,ypwrk
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
       common/gr/gr2
       common/dgr/dgr2,dgr,ddgr
 c
@@ -4049,14 +3962,14 @@ c local constants
 c arguments
       real(wp_), dimension(ndim) :: y,dery
 c local variables
-      integer iv,jv
+      integer :: iv,jv
       real(wp_) :: xx,yy,zz,yg2,anpl2,anpr2,del,dnl,duh,dan2sdnpl,
      . dan2sdxg,dan2sdyg,ddelnpl2,ddelnpl2x,ddelnpl2y,denom,derdel,
      . derdom,dfdiadnpl,dfdiadxg,dfdiadyg
       real(wp_), dimension(3) :: vgv,derdxv,danpldxv,derdnv,dbgr
 c common/external functions/variables
-      integer ierr
-      real(wp_) sox,gr2,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,anpl,
+      integer :: ierr
+      real(wp_) :: sox,gr2,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,anpl,
      . anpr,xg,yg,vgm,derdnm,dersdst
       real(wp_), dimension(3) :: dgr2,dgr,xv,anv,bv,derxg,deryg
       real(wp_), dimension(3,3) :: ddgr,derbv
@@ -4247,7 +4160,7 @@ c local variables
       real(wp_), dimension(3) :: dbtot,bvc
       real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv
 c common/external functions/variables
-      real(wp_) bres,brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr,
+      real(wp_) :: bres,brr,bphi,bzz,ajphi,btot,xg,yg,dxgdpsi,dpsidr,
      . dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv,psinv
       real(wp_), dimension(3) :: bv,derxg,deryg
       real(wp_), dimension(3,3) :: derbv
@@ -4768,7 +4681,7 @@ c arguments
 c local variables
       real(wp_) :: bzz,dbvcdc13,dbvcdc31
 c common/external functions/variables
-      real(wp_) dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
+      real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
 c
       common/derip1/dpsidr,dpsidz
       common/derip2/ddpsidrr,ddpsidzz,ddpsidrz
@@ -4819,7 +4732,7 @@ c
       use interp_eqprof, only : psia
       implicit none
 c common/external functions/variables
-      real(wp_) psinv,xg,xgcn,dxgdpsi,dens,ddenspsin
+      real(wp_) :: psinv,xg,xgcn,dxgdpsi,dens,ddenspsin
 c
       common/psival/psinv
       common/xgxg/xg
@@ -4978,9 +4891,11 @@ c
       use math, only : catand
       use graydata_flags, only : ipol
       use graydata_par, only : rwmax,psipol0,chipol0
+      use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
+     . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt
       implicit none
 c local constants
-      integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36
+      integer, parameter :: ndim=6,ndimm=3
 c local variables
       integer :: j,k,iproj,nfilp
       real(wp_) :: csth,snth,csps,snps,phiwrad,phirrad,csphiw,snphiw,
@@ -4995,32 +4910,15 @@ c local variables
       complex(wp_) :: sss,ddd,phic,qi1,qi2,tc,ts,qqxx,qqxy,qqyy,dqi1,
      . dqi2,dqqxx,dqqyy,dqqxy,d2qi1,d2qi2,d2qqxx,d2qqyy,d2qqxy
 c common/external functions/variables
-      integer :: nrayr,nrayth
       real(wp_) :: ak0,akinv,fhz,wcsi,weta,rcicsi,rcieta,phiw,phir,
      . anx0c,any0c,anz0c,x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi,
      . ddr11,psinv,dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi
-      real(wp_), dimension(jmx) :: dffiu,ddffiu
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: gri
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0
-      real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
-      complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt
 c
-      common/nray/nrayr,nrayth
       common/parwv/ak0,akinv,fhz
       common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
       common/anic/anx0c,any0c,anz0c
-      common/wrk/ywrk0,ypwrk0
-      common/grc/xc0,du10
-      common/fi/dffiu,ddffiu
       common/mirr/x00,y00,z00
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
-      common/evt/ext,eyt
       common/nplr/anpl,anpr
       common/psival/psinv
       common/parpl/brr,bphi,bzz,ajphi
@@ -5309,9 +5207,11 @@ c
      . cvdr=>degree,ui=>im
       use graydata_flags, only : ipol
       use graydata_par, only : rwmax,psipol0,chipol0
+      use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
+     . xc0=>xc,du10=>du1,dffiu,ddffiu,grad2,dgrad2v,gri,ggri,ext,eyt
       implicit none
 c local constants
-      integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36
+      integer, parameter :: ndim=6,ndimm=3
 c local variables
       integer :: j,k,iv,jv,iproj,nfilp
       real(wp_) :: csth,snth,csps,snps,phiwrad,csphiw,snphiw,dr,da,u,
@@ -5320,31 +5220,14 @@ c local variables
      . x0m,y0m,r0m,z0m,ancsi,aneta,ppcsi,ppeta,deltapol,qq,uu,vv
       real(wp_), dimension(ndim) :: ytmp,yptmp
 c common/external functions/variables
-      integer :: nrayr,nrayth
       real(wp_) :: wcsi,weta,rcicsi,rcieta,phiw,phir,anx0c,any0c,anz0c,
      . x00,y00,z00,dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,dens,ddens,
      . tekev,anpl,anpr,brr,bphi,bzz,ajphi
-      real(wp_), dimension(jmx) :: dffiu,ddffiu
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: gri
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0
-      real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
-      complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt
 c
-      common/nray/nrayr,nrayth
       common/parbeam/wcsi,weta,rcicsi,rcieta,phiw,phir
       common/anic/anx0c,any0c,anz0c
-      common/wrk/ywrk0,ypwrk0
-      common/grc/xc0,du10
-      common/fi/dffiu,ddffiu
       common/mirr/x00,y00,z00
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
-      common/evt/ext,eyt
       common/nplr/anpl,anpr
       common/psival/psinv
       common/parpl/brr,bphi,bzz,ajphi
@@ -5515,36 +5398,21 @@ c
       use const_and_precisions, only : wp_,izero,zero,one,pi,
      . cvdr=>degree
       use graydata_par, only : psipol0,chipol0
+      use beamdata, only : nrayr,nrayth,ywrk0=>ywrk,ypwrk0=>ypwrk,
+     . xc0=>xc,du10=>du1,grad2,dgrad2v,gri,ggri,yyrfl,ext,eyt
       implicit none
 c local constants
-      integer, parameter :: ndim=6,ndimm=3,jmx=31,kmx=36
+      integer, parameter :: ndim=6,ndimm=3
 c local variables
       integer :: j,k,iv,jv,iproj,nfilp
       real(wp_) :: x0,y0,z0,an20,an0,anx0,any0,anz0,vgradi,
      . r0,x0m,y0m,r0m,z0m,strfl11,qq,uu,vv
       real(wp_), dimension(ndim) :: ytmp,yptmp
 c common/external functions/variables
-      integer :: nrayr,nrayth
       real(wp_) :: dd,an2s,an2,fdia,bdotgr,ddi,ddr11,psinv,
      . dens,ddens,tekev,anpl,anpr,brr,bphi,bzz,ajphi
-      real(wp_), dimension(jmx,kmx) :: grad2
-      real(wp_), dimension(3,jmx,kmx) :: gri
-      real(wp_), dimension(ndim,jmx,kmx) :: ywrk0,ypwrk0
-      real(wp_), dimension(jmx,kmx,ndim) :: yyrfl
-      real(wp_), dimension(ndimm,jmx,kmx) :: xc0,du10,dgrad2v
-      real(wp_), dimension(3,3,jmx,kmx) :: ggri
-      complex(wp_), dimension(jmx,kmx,0:4) :: ext,eyt
 c
-      common/nray/nrayr,nrayth
-      common/wrk/ywrk0,ypwrk0
-      common/grc/xc0,du10
-      common/grad2jk/grad2
-      common/dgrad2jk/dgrad2v
-      common/gradjk/gri
-      common/ggradjk/ggri
       common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
-      common/yyrfl/yyrfl
-      common/evt/ext,eyt
       common/nplr/anpl,anpr
       common/psival/psinv
       common/parpl/brr,bphi,bzz,ajphi
@@ -5643,18 +5511,11 @@ c     ray power weigth coefficient q(j)
 c
       use const_and_precisions, only : wp_
       use graydata_par, only : rwmax
+      use beamdata, only : nrayr,nrayth,q
       implicit none
-c local constants
-      integer, parameter :: jmx=31
 c local variables
       integer :: j,k
       real(wp_) :: dr,r1,r2,summ,sm
-c common/external functions/variables
-      integer :: nrayr,nrayth
-      real(wp_), dimension(jmx) :: q
-c
-      common/qw/q
-      common/nray/nrayr,nrayth
 c
       dr=1.0_wp_
       if(nrayr.gt.1) dr=rwmax/dble(nrayr-1)
@@ -5754,9 +5615,9 @@ c
       use graydata_flags, only : dst
       use graydata_par, only : sgnbphi
       use graydata_anequil, only : rr0m
+      use beamdata, only : psjki,tauv,alphav,pdjki,ppabs,currj,didst,
+     . ccci,q,tau1v
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36,nmx=8000
 c arguments
       integer :: i,j,k
 c local variables
@@ -5766,18 +5627,7 @@ c local variables
 c common/external functions/variables
       real(wp_) :: p0mw,dersdst,alpha,effjcd,akim,tau0,psinv,
      . bmxi,bmni,fci
-      real(wp_), dimension(jmx) :: q
-      real(wp_), dimension(jmx,kmx) :: tau1v
-      real(wp_), dimension(jmx,kmx,nmx) :: psjki,tauv,alphav,pdjki,
-     . ppabs,currj,didst,ccci
 c
-      common/psjki/psjki
-      common/atjki/tauv,alphav
-      common/dpjjki/pdjki,currj
-      common/pcjki/ppabs,ccci
-      common/dcjki/didst
-      common/tau1v/tau1v
-      common/qw/q
       common/p0/p0mw
       common/dersdst/dersdst
       common/absor/alpha,effjcd,akim,tau0
@@ -5919,17 +5769,17 @@ c local constants
      . canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi)
       integer, parameter :: ksp=3
 c arguments
-      integer ieccd,nhmn,nhmx,ithn,iokhawa,ierr
+      integer :: ieccd,nhmn,nhmx,ithn,iokhawa,ierr
       real(wp_) :: yg,anpl,anprre,amu,Zeff,rbn,rbx,
      . fc,dens,psinv,effjcd
       complex(wp_) :: ex,ey,ez
 c local variables
-      integer nhn,njp,nlm,npar
+      integer :: nhn,njp,nlm,npar
       real(wp_) :: anum,denom,resp,resj,coullog,anucc,alams,
      . fp0s,pa,cst2
       real(wp_) :: chlm(nlmt)
       real(wp_), dimension(3+2*nlmt) :: eccdpar ! dimension(max(5,3+nlmt))
-c
+c common/external functions/variables
       real(wp_), external :: fjch,fjncl,fjch0
 c
       anum=0.0_wp_
@@ -6037,7 +5887,7 @@ c local variables
      . cst
       real(wp_), dimension(lw) :: w
       real(wp_), dimension(nfpp+necp) :: apar
-c
+c common/external functions/variables
       real(wp_), external :: fcur,fpp
 c
 c     EC power and current densities
@@ -6378,7 +6228,7 @@ c
       use const_and_precisions, only : wp_
       implicit none
 c arguments
-      real(wp_) upl,fjch0
+      real(wp_) :: upl,fjch0
       integer :: npar
       real(wp_), dimension(npar) :: extrapar
 c local variables
@@ -6453,9 +6303,9 @@ c arguments
       real(wp_) :: upl,fjncl
       real(wp_), dimension(npar) :: extrapar
 c local variables
+      integer :: nlm
       real(wp_) :: anpl,amu,ygn,zeff,fc,hb,gam,u2,u,upr2,
      . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp
-      integer :: nlm
 c
       anpl=extrapar(2)
       amu=extrapar(3)
@@ -6495,8 +6345,8 @@ c
 c local constants
       integer, parameter :: ksp=3
 c arguments
-      real(wp_) :: alam,fv,dfv
       integer :: nlmt
+      real(wp_) :: alam,fv,dfv
       real(wp_), dimension(nlmt) :: tlm,chlm
 c local variables
       integer :: nlm,ier
@@ -6519,22 +6369,18 @@ c
 c
       subroutine projxyzt(iproj,nfile)
       use const_and_precisions, only : wp_
+      use beamdata, only : ywrk,nrayr,nrayth
       implicit none
-c local constants
-      integer, parameter :: jmx=31,kmx=36
 c arguments
-      integer iproj,nfile
+      integer :: iproj,nfile
 c local variables
       integer :: jd,j,kkk,k
       real(wp_) :: dir,rtimn,rtimx,dx,dy,dz,dirx,diry,dirz,
      . csth1,snth1,csps1,snps1,xti,yti,zti,rti
 c common/external functions/variables
-      integer :: nrayr,nrayth,istep
+      integer :: istep
       real(wp_) :: psinv11,st
-      real(wp_), dimension(6,jmx,kmx) :: ywrk,ypwrk
 c
-      common/nray/nrayr,nrayth
-      common/wrk/ywrk,ypwrk
       common/psinv11/psinv11
       common/istep/istep
       common/ss/st
@@ -6612,9 +6458,10 @@ c
       use numint, only : simpson
       use graydata_flags, only : ipec,ieccd,nnd,dst
       use utils, only : locatex, locate, intlin
+      use beamdata, only : ppabs,ccci,psjki,nrayr,nrayth,nstep,iiv,pdjki
       implicit none
 c local constants
-      integer, parameter :: nndmx=5001,jmx=31,kmx=36,nmx=8000,llmx=21
+      integer, parameter :: nndmx=5001,llmx=21
       real(wp_), parameter :: rtbc=one
 c arguments
       real(wp_) :: pabs,currt
@@ -6631,24 +6478,17 @@ c local variables
      . rhotp,drhotp,stmx,pins_02,pins_05,pins_085,xrhot,currtka,rhop,
      . pinsr,ratjpli,ratjai,ratjbi,frhotor,frhopol,fdadrhot,
      . fdvdrhot,h,psin
-      real(wp_), dimension(nmx) :: xxi,ypt,yamp
+      real(wp_), dimension(nstep) :: xxi,ypt,yamp
       real(wp_), dimension(nndmx) :: ajphiv,dpdv,dvolt,darea,ratjav,
      . ratjbv,ratjplv,ajplv,ajcdav,ajcdbv,pins,currins,fi,rtab,rhotv
       real(wp_), dimension(0:nndmx) :: rtab1
 c common/external functions/variables
-      integer :: nrayr,nrayth,istep,index_rt
-      integer, dimension(jmx,kmx) :: iiv
-      real(wp_) :: alpha0,beta0,taumn,taumx,pabstot,currtot,
-     . psipol,chipol
-      real(wp_), dimension(jmx,kmx,nmx) :: psjki,ppabs,ccci,pdjki,currj
+      integer :: istep,index_rt
+      real(wp_) :: alpha0,beta0,taumn,taumx,pabstot,currtot,psipol,
+     . chipol
 c
-      common/nray/nrayr,nrayth
       common/istep/istep
       common/index_rt/index_rt
-      common/iiv/iiv
-      common/psjki/psjki
-      common/pcjki/ppabs,ccci
-      common/dpjjki/pdjki,currj
       common/angles/alpha0,beta0
       common/taumnx/taumn,taumx,pabstot,currtot
       common/polcof/psipol,chipol
@@ -6703,7 +6543,7 @@ c
         do k=1,kkk
           ise0=0
           ii=iiv(j,k)
-          if (ii.lt.nmx) then
+          if (ii.lt.nstep) then
             if(psjki(j,k,ii+1).ne.0.0_wp_) ii=ii+1
           end if
           idecr=-1
@@ -7079,7 +6919,7 @@ c arguments
 c local variables
       real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p
 c common/external functions/variables
-      real(wp_) anpl,anpr,xg,yg,sox
+      real(wp_) :: anpl,anpr,xg,yg,sox
 c
       common/nplr/anpl,anpr
       common/xgxg/xg