diff --git a/Makefile b/Makefile
index 131d69d..1e94bc7 100644
--- a/Makefile
+++ b/Makefile
@@ -1,61 +1,112 @@
-# Executable name
-EXE=gray
-LIBFILE=lib$(EXE).a
+# 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
 
-# Objects list
-OBJ=gray_main.o gray-externals.o grayl.o reflections.o green_func_p.o \
-    const_and_precisions.o itm_constants.o itm_types.o
+# library name
+# ------------
+LIBNAME=$(JLIBDIR)/libgray.a
 
-# Alternative search paths
-vpath %.f90 src
-vpath %.f src
+# List source and object files
+# ----------------------------
+FSRC=$(wildcard *.f *.f90)
+FOBJ0=$(FSRC:.f=.o)
+FOBJ=$(FOBJ0:.f90=.o)
+
+########################################
+# Set up compiler specific environment #
+########################################
 
-# Fortran compiler name and flags
-FC=gfortran
-FFLAGS=-O0 -g -fbounds-check -Wall
-#FFLAGS=-O3
+# set default compiler
+# --------------------
+FC=$(F90)
 
-DIRECTIVES = -DREVISION="'$(shell svnversion)'"
+# Set compiler debug and optimization flags
+# -----------------------------------------
+ifeq ("$(DBG)","")
+     FFLAGS= $(AUTO) -O3
+else
+     FFLAGS= $(AUTO) -O0 -g
+endif
+
+# Set include directories
+# -----------------------
+#INC=-I$(JCOMMON_STD)
+INC=
+
+#####################################
+# Set up RULES specific environment #
+#####################################
 
-all: library
+# suffixes
+# --------
+.SUFFIXES: .f90 .f .o .mod
 
-# Build executable from object files
-$(EXE): main.o $(OBJ)
-	$(FC) $(FFLAGS) -o $@ $^
+# 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 $(FFLAGS) $(INC)  $<
+
+# create rule for compiling f files
+# ---------------------------------
+%.o: %.f
+	$(FC) -c $(FFLAGS) $(INC)  $<
+
+# make 'realclean' option
+# -----------------------
+realclean:
+	rm -f *.o *.mod $(LIBNAME)
+
+# make 'clean' option
+# -------------------
+clean:
+	rm -f *.o
 
-$(LIBFILE): $(OBJ)
-	rm -f $@
-	ar -rv $@ $^
-
-# Dependencies on modules
-main.o: gray_main.o
+# Dependencies
+# ------------
 gray_main.o: gray-externals.o itm_types.o
 gray-externals.o: green_func_p.o reflections.o
 green_func_p.o: const_and_precisions.o
 const_and_precisions.o: itm_types.o itm_constants.o
 itm_constants.o: itm_types.o
 
-# General object compilation command
-%.o: %.f90
-	$(FC) $(FFLAGS) -c $<
-
+# Special rule to preprocess source file and include svn revision
+# ---------------------------------------------------------------
+DIRECTIVES = -DREVISION="'$(shell svnversion)'"
 gray-externals.o:gray-externals.f
-	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+	$(FC) -c -Mpreprocess $(DIRECTIVES) $(FFLAGS) $(INC)  $< 
 
-grayl.o:grayl.f
-	$(FC) $(FFLAGS) -c $<
-
-.PHONY: library clean install
-
-library: $(LIBFILE)
-
-# Remove output files
-clean:
-	rm -rf *.o *.mod $(EXE)
-
-install:
-	@if [ -f $(EXE) ]; then \
-	  cp $(EXE) ~/bin/; \
-	else \
-	  echo File $(EXE) does not exist. Run \'make\' first; \
-	fi
diff --git a/Makefile.standalone b/Makefile.standalone
new file mode 100644
index 0000000..131d69d
--- /dev/null
+++ b/Makefile.standalone
@@ -0,0 +1,61 @@
+# Executable name
+EXE=gray
+LIBFILE=lib$(EXE).a
+
+# Objects list
+OBJ=gray_main.o gray-externals.o grayl.o reflections.o green_func_p.o \
+    const_and_precisions.o itm_constants.o itm_types.o
+
+# Alternative search paths
+vpath %.f90 src
+vpath %.f src
+
+# Fortran compiler name and flags
+FC=gfortran
+FFLAGS=-O0 -g -fbounds-check -Wall
+#FFLAGS=-O3
+
+DIRECTIVES = -DREVISION="'$(shell svnversion)'"
+
+all: library
+
+# Build executable from object files
+$(EXE): main.o $(OBJ)
+	$(FC) $(FFLAGS) -o $@ $^
+
+$(LIBFILE): $(OBJ)
+	rm -f $@
+	ar -rv $@ $^
+
+# Dependencies on modules
+main.o: gray_main.o
+gray_main.o: gray-externals.o itm_types.o
+gray-externals.o: green_func_p.o reflections.o
+green_func_p.o: const_and_precisions.o
+const_and_precisions.o: itm_types.o itm_constants.o
+itm_constants.o: itm_types.o
+
+# General object compilation command
+%.o: %.f90
+	$(FC) $(FFLAGS) -c $<
+
+gray-externals.o:gray-externals.f
+	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+
+grayl.o:grayl.f
+	$(FC) $(FFLAGS) -c $<
+
+.PHONY: library clean install
+
+library: $(LIBFILE)
+
+# Remove output files
+clean:
+	rm -rf *.o *.mod $(EXE)
+
+install:
+	@if [ -f $(EXE) ]; then \
+	  cp $(EXE) ~/bin/; \
+	else \
+	  echo File $(EXE) does not exist. Run \'make\' first; \
+	fi
diff --git a/src/gray-externals.f b/src/gray-externals.f
index 6ac6d5e..8aa8a92 100644
--- a/src/gray-externals.f
+++ b/src/gray-externals.f
@@ -61,30 +61,30 @@ c     print all ray positions in local reference system
 c
       if(nrayr.gt.1) then
         iproj=1
-        nfilp=9
+        nfilp=609
         call projxyzt(iproj,nfilp)
       end if
 c
 c     print final results on screen
 c
-      print*,' '
-      print'(a,f9.4)','final step (s, ct, Sr)  = ',st
+      write(*,*)
+      write(*,'(a,f9.4)') 'final step (s, ct, Sr)  = ',st
       if(iwarm.gt.0) then
-      print '(a,2e12.5)','taumn, taumx  = ', taumn,taumx
+        write(*,'(a,2e12.5)') 'taumn, taumx  = ', taumn,taumx
       else
-      print '(a,2f9.4)','taumn, taumx  = ', zero,zero
+        write(*,'(a,2f9.4)') 'taumn, taumx  = ', zero,zero
       end if      
 c
-      print'(a,f9.4)','Pabs_tot (MW) = ',pabstot
+      write(*,'(a,f9.4)') 'Pabs_tot (MW) = ',pabstot
       currtka =currtot*1.0d3
-      print'(a,f9.4)','I_tot (kA)    = ',currtka
+      write(*,'(a,f9.4)') 'I_tot (kA)    = ',currtka
 c
       if (index_rt.eq.1) then
-      if(iequil.eq.2) write(6,*) 'EQUILIBRIUM CASE : ',filenmeqq
-      if(iequil.eq.1) write(6,*) 'ANALTYCAL EQUILIBRIUM'
-      if(iprof.eq.1) write(6,*)  'PROFILE file : ',filenmprf
-      if(iprof.eq.0) write(6,*)  'ANALTYCAL PROFILES'
-      if(ibeam.ge.1) write(6,*)  'LAUNCHER CASE : ',filenmbm
+      if(iequil.eq.2) write(*,*) 'EQUILIBRIUM CASE : ',filenmeqq
+      if(iequil.eq.1) write(*,*) 'ANALTYCAL EQUILIBRIUM'
+      if(iprof.eq.1) write(*,*)  'PROFILE file : ',filenmprf
+      if(iprof.eq.0) write(*,*)  'ANALTYCAL PROFILES'
+      if(ibeam.ge.1) write(*,*)  'LAUNCHER CASE : ',filenmbm
 
       end if
 c
@@ -161,7 +161,7 @@ c
           call gwork(j,k)
 c 
           if(ierr.gt.0) then
-            print*,'    IERR = ', ierr
+            write(*,*) '    IERR = ', ierr
             if(ierr.eq.97) then
 c              igrad=0
 c              ierr=0 
@@ -211,11 +211,11 @@ c         iov=2 first exit from plasma, iov=3 after 2nd entrance into plasma
                vva=sin(2.0d0*chipol0*cvdr)
                powa=0.5d0*(1.0d0+vva*vvin+uua*uuin+qqa*qqin)
 c              p0mw=p0mw*powa
-c              print*,' '
-c              print*,'Coupled power fraction =',powa
-c              print*,' '
-c              print*,'Input coupled power (MW) =',p0mw
-c              print*,' '
+c              write(*,*) ' '
+c              write(*,*) 'Coupled power fraction =',powa
+c              write(*,*) ' '
+c              write(*,*) 'Input coupled power (MW) =',p0mw
+c              write(*,*) ' '
             end if
             if (iop(j,k).eq.2.and.rrm.le.rcen.and.ipass.gt.1
      .          .and.ipolc.eq.1) then
@@ -229,7 +229,7 @@ c              print*,' '
               else
                 powrfl=0.5d0*(1.0d0+vvout*vvin2+uuout*uuin2+qqout*qqin2)
               end if
-              write(6,*) 'Reflected power fraction =',powrfl
+              write(*,*) 'Reflected power fraction =',powrfl
               iop(j,k)=3
               yyrfl(j,k,1)=xvrfl(1)
               yyrfl(j,k,2)=xvrfl(2)
@@ -276,10 +276,10 @@ c     print ray positions for j=nrayr in local reference system
 
       istpr=istpr+1
       if (istpr.eq.istpr0) then
-c        print*,'istep = ',i
+c        write(*,*) 'istep = ',i
         if(nrayr.gt.1) then
           iproj=0
-          nfilp=8
+          nfilp=608
           call projxyzt(iproj,nfilp)
         end if
         istpr=0
@@ -403,7 +403,7 @@ c     print dIds in A/m/W, Jphi and Jt in MA/m2, ki in m^-1
         tau11=taujk
         pt11=exp(-taujk)
 
-        write(4,99) stm,rrm,zzm,phideg,
+        write(604,99) stm,rrm,zzm,phideg,
      .  psi,rhot,dens11,tekev,
      .  bbr,bphi,bbz,ajphi*1.0d-6,sqrt(an2),anpl,akim*100,
      .  alpha11,tau11,pt11,ajcd11,dids11,
@@ -413,16 +413,16 @@ c        call polarcold(exf,eyif,ezf,elf,etf)
       end if
 c     central ray only end
 
-      if(k.eq.1.and.j.eq.nrayr) write(17,99) st,ddr11,ddr,ddi
+      if(k.eq.1.and.j.eq.nrayr) write(617,99) st,ddr11,ddr,ddi
 
 c     print dIds in MA/m, Jcd in MA/m^2, ppabs and pdjki in MW/m^3
       if(istpl.eq.istpl0)  then
         if(j.eq.nrayr) then
 c        if(j.eq.nrayr.and.tauv(j,k,i).le.taucr) then
-          write(33,111) i,j,k,stm,xxm,yym,rrm,zzm,
+          write(633,111) i,j,k,stm,xxm,yym,rrm,zzm,
      .          psjki(j,k,i),taujk,anpl,alphav(j,k,i),dble(index_rt)
         end if
-c      if(k.eq.nrayth) write(33,*) ' '
+c      if(k.eq.nrayth) write(633,*) ' '
       end if
 c
       return
@@ -439,27 +439,27 @@ c
 
       If(index_rt.eq.1) then
       
-      write(4,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '//
+      write(604,*)' #sst R z phi psi rhot ne Te Br Bphi Bz Jphi '//
      .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
-      write(8,*) ' #istep j k xt yt zt rt psin'
-      write(9,*) ' #istep j k xt yt zt rt psin'
-      write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
-      write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
-      write(12,*) ' #i sst psi w1 w2'
-      write(7,*)'#Icd Pa Jphip dPdVp '//
+      write(608,*) ' #istep j k xt yt zt rt psin'
+      write(609,*) ' #istep j k xt yt zt rt psin'
+      write(617,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
+      write(633,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
+      write(612,*) ' #i sst psi w1 w2'
+      write(607,*)'#Icd Pa Jphip dPdVp '//
      .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
      .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
-      write(48,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
+      write(648,*) '#psi rhot Jphi Jcdb dPdV Icdins Pins P% index_rt'
 
       else
 
 c      If(index_rt.eq.3) then
-      write(4,*) ' '
-      write(8,*) ' '
-      write(9,*) ' '
-      write(17,*) ' '
-      write(12,*) ' '
-      write(48,*) ' '
+      write(604,*) ' '
+      write(608,*) ' '
+      write(609,*) ' '
+      write(617,*) ' '
+      write(612,*) ' '
+      write(648,*) ' '
 
 c      end if
       end if
@@ -541,15 +541,15 @@ c
       common/angles/alpha0,beta0
       common/scal/iscal
 c
-      open(2,file='gray.data',status= 'unknown')
+      open(602,file='gray.data',status= 'unknown')
 c
 c     alpha0, beta0 (cartesian) launching angles
 c     fghz      wave frequency (GHz)
 c     p0mw      injected power (MW)
 c
-      read(2,*) alpha0,beta0
-      read(2,*) fghz
-      read(2,*) p0mw
+      read(602,*) alpha0,beta0
+      read(602,*) fghz
+      read(602,*) p0mw
 c     power value overwritten with value passed as argument (in W)
       p0mw=powin*1.d-6
 c
@@ -558,36 +558,36 @@ c     nrayth      number of rays in angular direction
 c     rwmax      normalized maximum radius of beam power
 c     rwmax=1 -> last ray at radius = waist
 c
-      read(2,*) nrayr,nrayth,rwmax
+      read(602,*) nrayr,nrayth,rwmax
       if(nrayr.eq.1) nrayth=1
 c
 c     x00,y00,z00 coordinates of launching point
 c
-      read(2,*) x00,y00,z00
+      read(602,*) x00,y00,z00
 c
 c     beams parameters in local reference system
 c     w0 -> waist, d0 -> waist distance from launching point
 c     phiw  angle of beam ellipse
 c
-      read(2,*) w0csi,w0eta,d0csi,d0eta,phiw
+      read(602,*) w0csi,w0eta,d0csi,d0eta,phiw
 c
 c     ibeam=0  :read data for beam as above
 c     ibeam=1  :read data from file simple astigmatic beam
 c     ibeam=2  :read data from file general astigmatic beam
-      read(2,*) ibeam
-      read(2,*) filenmbm
+      read(602,*) ibeam
+      read(602,*) filenmbm
 c
 c     iequil=0  :vacuum
 c     iequil=1  :analytical equilibrium
 c     iequil=2  :read eqdsk
 c     ixp=0,-1,+1  : no X point , bottom/up X point
 c
-      read(2,*) iequil,ixp
+      read(602,*) iequil,ixp
 c
 c     iprof=0  :analytical density and temp. profiles
 c     iprof>0  :numerical density and temp. profiles
 c
-      read(2,*) iprof
+      read(602,*) iprof
 c
 c     iwarm=0  :no absorption and cd
 c     iwarm=1  :weakly relativistic absorption
@@ -595,18 +595,18 @@ c     iwarm=2  :relativistic absorption, n<1 asymptotic expansion
 c     iwarm=3  :relativistic absorption, numerical integration
 c     ilarm    :order of larmor expansion
 c
-      read(2,*) iwarm,ilarm
+      read(602,*) iwarm,ilarm
 c      if(iwarm.gt.2) iwarm=2
 c
 c     ieccd 0/1 NO/YES  ECCD calculation ieccd>0 different CD models
 c
-      read(2,*) ieccd
+      read(602,*) ieccd
       if (ieccd.ge.10)     call Setup_SpitzFunc
 c
 c     ipec=0/1  :pec profiles grid in psi/rhop
 c     nnd       :number of grid steps for pec profiles +1
 c
-      read(2,*) ipec,nnd
+      read(602,*) ipec,nnd
 c
 c     dst       integration step
 c     nstep      maximum number of integration steps
@@ -620,15 +620,15 @@ c     ipass=1/2 1 or 2 passes into plasma
 c     iox=1/2    OM/XM
 c     idst=0/1/2  0 integration in s, 1 integr. in ct, 2 integr. in Sr
 c
-      read(2,*) dst,nstep,istpr0,istpl0,idst
-      read(2,*) igrad,ipass,rwallm
-      read(2,*) iox,psipol0,chipol0
+      read(602,*) dst,nstep,istpr0,istpl0,idst
+      read(602,*) igrad,ipass,rwallm
+      read(602,*) iox,psipol0,chipol0
 c
 c     ipsinorm  0 standard EQDSK format, 1 format Portone summer 2004
 c     sspl      spline parameter for psi interpolation
 c
-      read(2,*) filenmeqq
-      read(2,*) ipsinorm,sspl
+      read(602,*) filenmeqq
+      read(602,*) ipsinorm,sspl
 c     psin grid is normalized
       ipsinorm=1
 c
@@ -637,18 +637,18 @@ c     scaling adopted: beta=const, qpsi=const, nustar=const
 c     iscal    ne Te scaling 0: nustar=const,  1: n_greenw=const; 2 no rescaling
 c     factT factn   factor for Te&ne scaling 
 c
-      read(2,*) factb, iscal,factt,factn
+      read(602,*) factb, iscal,factt,factn
       if (factb.eq.0.0d0) factb=1.0d0
 c
 c     psbnd     value of psi ( > 1 ) of density boundary
 c
-      read(2,*) filenmprf
-      read(2,*) psdbnd
+      read(602,*) filenmprf
+      read(602,*) psdbnd
       if(iprof.eq.0) psdbnd=1.0d0
 c
 c     signum of toroidal B and I
 c     icocos index for equilibrium from COCOS - O. Sauter Feb 2012
-      read(2,*) sgnbphi,sgniphi,icocos
+      read(602,*) sgnbphi,sgniphi,icocos
 c     data are already in the coordinate system used by gray
       icocos=3
 c
@@ -656,37 +656,37 @@ c     read parameters for analytical density and temperature
 c     profiles if iprof=0
 c
       if(iprof.eq.0) then
-        read(2,*) dens0,aln1,aln2
-        read(2,*) te0,dte0,alt1,alt2
+        read(602,*) dens0,aln1,aln2
+        read(602,*) te0,dte0,alt1,alt2
       else
-        read(2,*) dummy,dummy,dummy
-        read(2,*) dummy,dummy,dummy,dummy
+        read(602,*) dummy,dummy,dummy
+        read(602,*) dummy,dummy,dummy,dummy
       end if
-      read(2,*) zeffan
+      read(602,*) zeffan
 c
 c     read parameters for analytical simple cylindical equilibrium
 c     if iequil=0,1
 
       if(iequil.lt.2) then
-        read(2,*) rr0,zr0,rpa
-        read(2,*) b0
-        read(2,*) q0,qa,alq
+        read(602,*) rr0,zr0,rpa
+        read(602,*) b0
+        read(602,*) q0,qa,alq
         rr0m=rr0/1.0d2
         zr0m=zr0/1.0d2
         rpam=rpa/1.0d2
       else
-        read(2,*) dummy,dummy,dummy
-        read(2,*) dummy
-        read(2,*) dummy,dummy,dummy
+        read(602,*) dummy,dummy,dummy
+        read(602,*) dummy
+        read(602,*) dummy,dummy,dummy
       end if
 c
-      close(unit=2)
+      close(unit=602)
 c
       if(nrayr.eq.1) igrad=0
       if (nrayr.lt.5) then
         igrad=0
-        print*,'    nrayr < 5 ! => OPTICAL CASE ONLY'
-        print*,' '
+        write(*,*) '    nrayr < 5 ! => OPTICAL CASE ONLY'
+        write(*,*)
       end if
 c
       fhz=fghz*1.0d9
@@ -730,13 +730,13 @@ c
         end if
       end if
 c
-      print'(a,2f8.3)','alpha0, beta0 = ',alpha0,beta0
+      write(*,'(a,2f8.3)') 'alpha0, beta0 = ',alpha0,beta0
 c
       r00=sqrt(x00**2+y00**2)
       phi0=acos(x00/r00)
       if(y00.lt.0) phi0=-phi0
-      print'(a,4f8.3)','x00, y00, R00, z00 = ',x00,y00,r00,z00
-      print*,' '
+      write(*,'(a,4f8.3)') 'x00, y00, R00, z00 = ',x00,y00,r00,z00
+      write(*,*)
 c
 c      anx0c=(anr0c*x00-anphi0c*y00)/r00
 c      any0c=(anr0c*y00+anphi0c*x00)/r00
@@ -806,7 +806,7 @@ c     set simple limiter as two cylindrical walls at rwallm and r00
         ipass=abs(ipass)
       end if
 
-      nfil=78
+      nfil=638
       
       open(nfil,file='headers.txt',status= 'unknown')
       call date_and_time(wdat,wtim)
@@ -872,7 +872,7 @@ c
 c
 c     print circular magnetic surfaces iequil=1
 c
-      write(71,*) '#i psin R z'
+      write(631,*) '#i psin R z'
       dal=2.0d-2*pi
       drn=0.1d0
       do i=1,10
@@ -882,21 +882,21 @@ c
           dzzm=rpam*rni*sin((k-1)*dal)
           rrm=rr0m+drrm
           zzm=zr0m+dzzm
-          write(71,111) i,rni,rrm,zzm
+          write(631,111) i,rni,rrm,zzm
         end do
-        write(71,*) ' '
-        write(71,*) ' '
+        write(631,*) ' '
+        write(631,*) ' '
       end do
 c
 c     print resonant B iequil=1
 c
-      write(70,*)'#i Btot  R  z'
+      write(630,*)'#i Btot  R  z'
       rres=b0*rr0m/bres
       zmx=zr0m+rpam
       zmn=zr0m-rpam
       do i=1,21
           zzres=zmn+(i-1)*(zmx-zmn)/2.0d1
-          write(70,111) i,bres,rres,zzres
+          write(630,111) i,bres,rres,zzres
       end do
 
       return
@@ -933,7 +933,7 @@ c     ibeam=2 general astigmatic beam
 c
 c     initial beam data are measured in mm -> transformed to cm
 c
-      nfbeam=97
+      nfbeam=603
       lbm=length(filenmbm)
       open(file=filenmbm(1:lbm)//'.txt',status= 'unknown',unit=nfbeam)
 c
@@ -968,32 +968,32 @@ c
       end do
 c
       iopt=0
-      call difcs(alphastv,betastv,nisteer,iopt,cbeta,ier)
-      call difcs(alphastv,waist1v,nisteer,iopt,cwaist1,ier)
-      call difcs(alphastv,rci1v,nisteer,iopt,crci1,ier)
-      call difcs(alphastv,waist2v,nisteer,iopt,cwaist2,ier)
-      call difcs(alphastv,rci2v,nisteer,iopt,crci2,ier)
-      call difcs(alphastv,phi1v,nisteer,iopt,cphi1,ier)
-      call difcs(alphastv,phi2v,nisteer,iopt,cphi2,ier)
-      call difcs(alphastv,x00v,nisteer,iopt,cx0,ier)
-      call difcs(alphastv,y00v,nisteer,iopt,cy0,ier)
-      call difcs(alphastv,z00v,nisteer,iopt,cz0,ier)
+      call difcsg(alphastv,betastv,nisteer,iopt,cbeta,ier)
+      call difcsg(alphastv,waist1v,nisteer,iopt,cwaist1,ier)
+      call difcsg(alphastv,rci1v,nisteer,iopt,crci1,ier)
+      call difcsg(alphastv,waist2v,nisteer,iopt,cwaist2,ier)
+      call difcsg(alphastv,rci2v,nisteer,iopt,crci2,ier)
+      call difcsg(alphastv,phi1v,nisteer,iopt,cphi1,ier)
+      call difcsg(alphastv,phi2v,nisteer,iopt,cphi2,ier)
+      call difcsg(alphastv,x00v,nisteer,iopt,cx0,ier)
+      call difcsg(alphastv,y00v,nisteer,iopt,cy0,ier)
+      call difcsg(alphastv,z00v,nisteer,iopt,cz0,ier)
 c
       if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
-        call locate(alphastv,nisteer,alpha0,k)
+        call vlocate(alphastv,nisteer,alpha0,k)
         dal=alpha0-alphastv(k)
-        betst=spli(cbeta,nisteer,k,dal)
-        x00=spli(cx0,nisteer,k,dal)
-        y00=spli(cy0,nisteer,k,dal)
-        z00=spli(cz0,nisteer,k,dal)
-        wcsi=spli(cwaist1,nisteer,k,dal)
-        weta=spli(cwaist2,nisteer,k,dal)
-        rcicsi=spli(crci1,nisteer,k,dal)
-        rcieta=spli(crci2,nisteer,k,dal)
-        phiw=spli(cphi1,nisteer,k,dal)
-        phir=spli(cphi2,nisteer,k,dal)
+        betst=fspli(cbeta,nisteer,k,dal)
+        x00=fspli(cx0,nisteer,k,dal)
+        y00=fspli(cy0,nisteer,k,dal)
+        z00=fspli(cz0,nisteer,k,dal)
+        wcsi=fspli(cwaist1,nisteer,k,dal)
+        weta=fspli(cwaist2,nisteer,k,dal)
+        rcicsi=fspli(crci1,nisteer,k,dal)
+        rcieta=fspli(crci2,nisteer,k,dal)
+        phiw=fspli(cphi1,nisteer,k,dal)
+        phir=fspli(cphi2,nisteer,k,dal)
       else
-        print*,'    alpha0 outside table range !!!'
+        write(*,*) '    alpha0 outside table range !!!'
         if(alpha0.ge.alphastv(nisteer)) ii=nisteer
         if(alpha0.le.alphastv(1)) ii=1
         betst=betastv(ii)
@@ -1125,9 +1125,9 @@ c     psi function
 c
 c      do j=1,nz
 c        do i=1,nr
-c          write(80,2021) rv(i),zv(j),psi(i,j)
+c          write(620,2021) rv(i),zv(j),psi(i,j)
 c        enddo
-c        write(80,*) ' '
+c        write(620,*) ' '
 c      enddo
 
       do j=1,nz
@@ -1162,9 +1162,9 @@ c
             ij=nz*(i-1)+j
             psin(i,j)=ffvpsi(ij)
             psi(i,j)=psin(i,j)*psia
-c            write(79,2021) rv(i),zv(j),psin(i,j)
+c            write(619,2021) rv(i),zv(j),psin(i,j)
           enddo
-c            write(79,*) ' '
+c            write(619,*) ' '
         enddo
       end if
 2021  format(5(1x,e16.9))
@@ -1254,7 +1254,7 @@ c
       call points_ox(rmaxis,zmaxis,rmop,zmop,psinoptmp,info)
       rmaxis=rmop
       zmaxis=zmop
-      print'(a,2f8.4,es12.5)','O-point',rmop,zmop,psinoptmp
+      write(*,'(a,2f8.4,es12.5)') 'O-point',rmop,zmop,psinoptmp
 c
 c     search for X-point   if ixp not = 0
 c
@@ -1264,7 +1264,7 @@ c
           z10=zbmin
           call points_ox(r10,z10,rxp,zxp,psinxptmp,info)
           if(psinxp.ne.-1.0d0) then
-            print'(a,2f8.4,es12.5)','X-point',rxp,zxp,psinxptmp
+            write(*,'(a,2f8.4,es12.5)') 'X-point',rxp,zxp,psinxptmp
             rbmin=rxp
             zbmin=zxp
             psinop=psinoptmp
@@ -1278,14 +1278,14 @@ c
             zbmax=z1
           else
             ixp=0
-c            print'(a)','no X-point'
+c            write(*,'(a)') 'no X-point'
            end if
         else
           r10=rmop
           z10=zbmax
           call points_ox(r10,z10,rxp,zxp,psinxptmp,info)
           if(psinxp.ne.-1.0d0) then
-            print'(a,2f8.4,e16.8)','X-point',rxp,zxp,psinxptmp
+            write(*,'(a,2f8.4,e16.8)') 'X-point',rxp,zxp,psinxptmp
             zbmax=zxp
             psinop=psinoptmp
             psinxp=psinxptmp
@@ -1296,7 +1296,7 @@ c            print'(a)','no X-point'
             zbmin=z1
           else
             ixp=0
-c           print'(a)','no X-point'
+c           write(*,'(a)') 'no X-point'
            end if
         end if
       end if
@@ -1314,24 +1314,24 @@ c
         call points_tgo(r10,z10,r1,z1,psin1,info)
         zbmin=z1
         rbmin=r1
-        print'(a,4f8.4)','no X-point  ',rbmin,zbmin,rbmax,zbmax
+        write(*,'(a,4f8.4)') 'no X-point  ',rbmin,zbmin,rbmax,zbmax
       end if
-      print*,' '
+      write(*,*) ' '
 
 c     compute B_toroidal on axis
 
       btaxis=fpol(1)/rmaxis
       btrcen=abs(btrcen)*factb
-      print'(a,f8.4)','factb =   ',factb
-      print'(a,f8.4)','|BT_centr|= ',btrcen
-      print'(a,f8.4)','|BT_axis| = ',abs(btaxis)
+      write(*,'(a,f8.4)') 'factb =   ',factb
+      write(*,'(a,f8.4)') '|BT_centr|= ',btrcen
+      write(*,'(a,f8.4)') '|BT_axis| = ',abs(btaxis)
 
 c     compute normalized rho_tor from eqdsk q profile
       call rhotor(nrho)
       phitedge=abs(psia)*rhotsx*2*pi
       rrtor=sqrt(phitedge/abs(btrcen)/pi)
       call rhopol
-c      print*,rhotsx,phitedge,rrtor,abs(psia)
+c      write(*,*) rhotsx,phitedge,rrtor,abs(psia)
 
 c     compute flux surface averaged quantities
 
@@ -1343,7 +1343,7 @@ c     compute flux surface averaged quantities
 c      ipr=1
 c      call contours_psi(1.0d0,np,rcon,zcon,ipr)
 c      do ii=1,2*np+1
-c        write(52,*) rcon(ii), zcon(ii)
+c        write(632,*) rcon(ii), zcon(ii)
 c      end do
 c
 
@@ -1359,13 +1359,13 @@ c     locate  psi surface for q=1.5 and q=2
       if (q15.gt.qmin.and.q15.lt.qmax) then
         call surfq(q15,psi15)
         rhot15=frhotor(psi15)
-        print'(3(a,f8.5))','psi_1.5 = ',psi15,'  rhop_1.5 = '
+        write(*,'(3(a,f8.5))') 'psi_1.5 = ',psi15,'  rhop_1.5 = '
      .  ,sqrt(psi15),'  rhot_1.5 = ',rhot15
       end if
       if (q2.gt.qmin.and.q2.lt.qmax) then
         call surfq(q2,psi2)
         rhot2=frhotor(psi2)
-        print'(3(a,f8.5))','psi_2   = ',psi2,'  rhop_2   = '
+        write(*,'(3(a,f8.5))') 'psi_2   = ',psi2,'  rhop_2   = '
      .   ,sqrt(psi2),'  rhot_2   = ',rhot2
       end if
 c
@@ -1389,7 +1389,7 @@ c
       tol = sqrt(dpmpar(1))
       call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
       if(info.gt.1)  then
-        print'(a,i2,a,2f8.4)','    info subr points_ox =',info,
+        write(*,'(a,i2,a,2f8.4)') '    info subr points_ox =',info,
      .   '  O/X coord.',xvec
       end if
       rf=xvec(1)
@@ -1472,7 +1472,7 @@ c
       common/eqnn/nr,nz,nrho,npp,nintp
       common/dens/dens,ddens
 c
-      write(55,*) ' #psi rhot ne Te q Jphi'
+      write(645,*) ' #psi rhot ne Te q Jphi'
       psin=0.0d0
       rhop=0.0d0
       rhot=0.0d0
@@ -1481,7 +1481,7 @@ c
       te=temp(psin)
       qq=qpsi(1)
 c
-      write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
+      write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6
 c
       nst=nrho
       do i=2,nst
@@ -1491,13 +1491,13 @@ c
         call density(psin)
         te=temp(psin)
 c
-        call locate(psinr,nrho,psin,ips)
+        call vlocate(psinr,nrho,psin,ips)
         ips=min(max(ips,1),nrho-1)
         call intlin(psinr(ips),qpsi(ips),psinr(ips+1),qpsi(ips+1),
      .            psin,qq)
         rhot=frhotor(psin)
         call tor_curr_psi(psin,ajphi)
-        write(55,111) psin,rhot,dens,te,qq,ajphi*1.d-6
+        write(645,111) psin,rhot,dens,te,qq,ajphi*1.d-6
       end do
 c
       return
@@ -1519,7 +1519,7 @@ c
 c
 c     locate  psi surface for q=qval
 c
-      call locate(qpsi,nrho,qval,i1)
+      call vlocate(qpsi,nrho,qval,i1)
       call intlin(qpsi(i1),psinr(i1),qpsi(i1+1),psinr(i1+1),
      .  qval,psival)
       ipr=1
@@ -1554,23 +1554,23 @@ c
           btotal(j,k)=sqrt(bbr**2+bbz**2+bbphi**2)
           if(btotal(j,k).ge.btmx) btmx=btotal(j,k)
           if(btotal(j,k).le.btmn) btmn=btotal(j,k)
-c          write(90,113) j,rrj,zzk,btotal(j,k)
+c          write(640,113) j,rrj,zzk,btotal(j,k)
         enddo
-c          write(90,*) ' '
+c          write(640,*) ' '
        enddo
 c
 c     compute Btot=Bres and Btot=Bres/2
 c
-      write(70,*)'#i Btot  R  z'
+      write(630,*)'#i Btot  R  z'
       do n=1,5
         bbb=bres/dble(n)
         if (bbb.ge.btmn.and.bbb.le.btmx) then
           call cniteq(bbb,nconts,ncpts,nctot,rrcb,zzcb,1)
           do inc=1,nctot
-            write(70,113) inc,bbb,rrcb(inc),zzcb(inc)
+            write(630,113) inc,bbb,rrcb(inc),zzcb(inc)
           end do
         end if
-        write(70,*) ' '
+        write(630,*) ' '
       end do
 c
       return
@@ -1650,12 +1650,12 @@ c
 c
       call splev(tfn,nsfd,cfn,3,psrad,densi,npp,ier)
       nu=1
-      call splder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
+      call gsplder(tfn,nsfd,cfn,3,nu,psrad,ddensi,npp,wrkfd,ier)
       dnpp=densi(npp)
       ddnpp=ddensi(npp)
 c
       nu=2
-      call splder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
+      call gsplder(tfn,nsfd,cfn,3,nu,psrad,d2densi,npp,wrkfd,ier)
       d2dnpp=d2densi(npp)
 
         if(derad(npp).eq.0.0d0) then
@@ -1711,7 +1711,7 @@ c
 c     spline interpolation of rhotor
 c
       iopt=0
-      call difcs(psinr,rhotnr,nr,iopt,crhot,ier)
+      call difcsg(psinr,rhotnr,nr,iopt,crhot,ier)
       return
       end
 
@@ -1726,7 +1726,7 @@ c
       if(irt.eq.0) irt=1
       if(irt.eq.nrho) irt=nrho-1
       dps=psi-psinr(irt)
-      fq_eq=spli(cq,nrho,irt,dps)
+      fq_eq=fspli(cq,nrho,irt,dps)
       return
       end
 
@@ -1742,7 +1742,7 @@ c
       if(irt.eq.0) irt=1
       if(irt.eq.nrho) irt=nrho-1
       dps=psi-psinr(irt)
-      frhotor_eq=spli(crhot,nrho,irt,dps)
+      frhotor_eq=fspli(crhot,nrho,irt,dps)
       return
       end
 
@@ -1765,7 +1765,7 @@ c
       if(ip.eq.0) ip=1
       if(ip.eq.nintp) ip=nintp-1
       dps=rpsi-rpstab(ip)
-      frhotor_av=spli(crhotq,nintp,ip,dps)
+      frhotor_av=fspli(crhotq,nintp,ip,dps)
       return
       end
       
@@ -1798,10 +1798,10 @@ c     spline interpolation of rhopol versus rhotor
       kspl=3
       call curfit(iopt,nnr,rhot,rhop,wp,xb,xe,kspl,ss,nrest,nsrp,
      . trp,crp,rp,wrkp,lwrkp,iwrkp,ier)
-c      print*,ier
+c      write(*,*) ier
            call splev(trp,nsrp,crp,3,rhot,rhopi,nnr,ier)
       do i=1,nnr
-        write(54,*) rhop(i),rhot(i),rhopi(i)
+        write(644,*) rhop(i),rhot(i),rhopi(i)
       end do 
            
       return
@@ -2046,8 +2046,8 @@ c
       do ic=2,np
         zc=zlw+(zup-zlw)*(1.0d0-cos(th*(ic-1)))/2.0d0
         iopt=1
-        call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
-        if(ier.gt.0) print*,'    profil =',ier
+        call sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
+        if(ier.gt.0) write(*,*) '    sprofil =',ier
         val=h*psiant+psinop
         call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)      
         if (zeroc(1).gt.rwallm) then
@@ -2064,10 +2064,10 @@ c
       end do
       if (ipr.gt.0) then
       do ii=1,2*np+1
-        write(71,111) ii,h,rcn(ii),zcn(ii)
+        write(631,111) ii,h,rcn(ii),zcn(ii)
       end do
-        write(71,*)
-        write(71,*)
+        write(631,*)
+        write(631,*)
       end if
       return
 111   format(i6,12(1x,e12.5))
@@ -2133,7 +2133,7 @@ c
 
 c     computation of flux surface averaged quantities
 
-      write(71,*)' #i psin R  z'
+      write(631,*)' #i psin R  z'
 
       nintp=nnintp
       ninpr=(nintp-1)/10
@@ -2313,7 +2313,7 @@ c     ratio_pltor  = Jcd_||/J_phi         Jcd_|| = <Jcd.b>
         dvdrhotv(jp)=phitedge*frhotor_eq(height2)/fq_eq(height2)
      .   *dvdpsi/pi
 c        
-c      write(57,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)
+c      write(647,99) sqrt(pstab(jp)),pstab(jp),riav,dvdpsi,area,vvol(jp)
 
 c       computation of rhot from calculated q profile 
 
@@ -2365,7 +2365,7 @@ c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
         end do
       end do
 
-      write(56,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
+      write(646,*)' #psi rhot_eq rhot_av |<B>| |Bmx| |Bmn|'//
      .'  Area Vol |I_pl| <J_phi> qq fc ratioJa ratioJb'
 
       qqv(1)=qqv(2)
@@ -2378,7 +2378,7 @@ c  ffhlam = Bmn/Bmx/fc integral_lambda^1 dlam/<sqrt(1-lam*B(rhop)/Bmx)>
               pstab(jp)=1.0d0
         end if
         rhot_eq=frhotor_eq(pstab(jp))
-        write(56,99) pstab(jp),rhot_eq,rhotqv(jp),
+        write(646,99) pstab(jp),rhot_eq,rhotqv(jp),
      .    bav(jp),bmxpsi(jp),bmnpsi(jp),varea(jp),vvol(jp), 
      .    vcurrp(jp),vajphiav(jp),qqv(jp),ffc(jp)
      .    ,vratja(jp),vratjb(jp)
@@ -2391,31 +2391,31 @@ c     used for computations of dP/dV and J_cd
 c     spline coefficients of rhot
 
         iopt=0
-        call difcs(rpstab,vvol,nintp,iopt,cvol,ier)
+        call difcsg(rpstab,vvol,nintp,iopt,cvol,ier)
         iopt=0
-        call difcs(rpstab,rbav,nintp,iopt,crbav,ier)
+        call difcsg(rpstab,rbav,nintp,iopt,crbav,ier)
         iopt=0
-        call difcs(rpstab,rri,nintp,iopt,crri,ier)
+        call difcsg(rpstab,rri,nintp,iopt,crri,ier)
         iopt=0
-        call difcs(rpstab,bmxpsi,nintp,iopt,cbmx,ier)
+        call difcsg(rpstab,bmxpsi,nintp,iopt,cbmx,ier)
         iopt=0
-        call difcs(rpstab,bmnpsi,nintp,iopt,cbmn,ier)
+        call difcsg(rpstab,bmnpsi,nintp,iopt,cbmn,ier)
         iopt=0
-        call difcs(rpstab,vratja,nintp,iopt,cratja,ier)
+        call difcsg(rpstab,vratja,nintp,iopt,cratja,ier)
         iopt=0
-        call difcs(rpstab,vratjb,nintp,iopt,cratjb,ier)
+        call difcsg(rpstab,vratjb,nintp,iopt,cratjb,ier)
         iopt=0
-        call difcs(rpstab,vratjpl,nintp,iopt,cratjpl,ier)
+        call difcsg(rpstab,vratjpl,nintp,iopt,cratjpl,ier)
         iopt=0
-        call difcs(rpstab,varea,nintp,iopt,carea,ier)
+        call difcsg(rpstab,varea,nintp,iopt,carea,ier)
         iopt=0
-        call difcs(rpstab,ffc,nintp,iopt,cfc,ier)
+        call difcsg(rpstab,ffc,nintp,iopt,cfc,ier)
         iopt=0
-        call difcs(rpstab,rhotqv,nintp,iopt,crhotq,ier)
+        call difcsg(rpstab,rhotqv,nintp,iopt,crhotq,ier)
         iopt=0
-        call difcs(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier)
+        call difcsg(rpstab,dadrhotv,nintp,iopt,cdadrhot,ier)
         iopt=0
-        call difcs(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier)
+        call difcsg(rpstab,dvdrhotv,nintp,iopt,cdvdrhot,ier)
 
 c   spline interpolation of H(lambda,rhop) and dH/dlambda
 
@@ -2446,7 +2446,7 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       if(ip.eq.0) ip=1
       if(ip.eq.nintp) ip=nintp-1
       dps=rpsi-rpstab(ip)
-      fdadrhot=spli(cdadrhot,nintp,ip,dps)
+      fdadrhot=fspli(cdadrhot,nintp,ip,dps)
       return
       end
 
@@ -2461,7 +2461,7 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       if(ip.eq.0) ip=1
       if(ip.eq.nintp) ip=nintp-1
       dps=rpsi-rpstab(ip)
-      fdvdrhot=spli(cdvdrhot,nintp,ip,dps)
+      fdvdrhot=fspli(cdvdrhot,nintp,ip,dps)
       return
       end
 
@@ -3435,7 +3435,7 @@ c
      .     rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
       psinv=(ffspl(1)-psinop)/psiant
       if(psinv.lt.0.0d0)
-     .   print'(a,3e12.4)', '    psin < 0 , R, z  ',psinv,rpsim,zpsim
+     .   write(*,'(a,3e12.4)') '    psin < 0 , R, z  ',psinv,rpsim,zpsim
 c
       nur=1
       nuz=0
@@ -3492,7 +3492,7 @@ c
         call splev(tfp,nsft,cfp,3,rrs,ffspl,1,ier)
         fpolv=ffspl(1)
         nu=1
-        call splder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
+        call gsplder(tfp,nsft,cfp,3,nu,rrs,ffspl,1,wrkfd,ier)
         dfpolv=ffspl(1)
         ffpv=fpolv*dfpolv/psia
       end if
@@ -3550,8 +3550,8 @@ c
 c
       iopt=1
       zc=zmaxis
-      call profil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
-      if(ier.gt.0) print*,'    profil =',ier
+      call sprofil(iopt,tr,nsr,tz,nsz,cc,kspl,kspl,zc,nsr,czc,ier)
+      if(ier.gt.0) write(*,*) '    sprofil =',ier
       val=h*psiant+psinop
       call sproota(val,tr,nsr,czc,zeroc,mest,m,ier)
       r1=zeroc(1)
@@ -3628,12 +3628,12 @@ c
           dens=ffs(1)
           nu=1
           ier=0
-          call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
+          call gsplder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
           ddens=ffs(1)
-          if(ier.gt.0) print*,ier
+          if(ier.gt.0) write(*,*) ier
           if(abs(dens).lt.1.0d-10) dens=0.0d0
         end if
-        if(dens.lt.0.0d0) print*,'    DENSITY NEGATIVE',dens
+        if(dens.lt.0.0d0) write(*,*) '    DENSITY NEGATIVE',dens
 c
       end if
       return
@@ -3658,11 +3658,11 @@ c
         proft=(1.0d0-arg**alt1)**alt2
         temp=(te0-dte0)*proft+dte0
       else
-        call locate(psrad,npp,arg,k)
+        call vlocate(psrad,npp,arg,k)
         if(k.eq.0) k=1
         if(k.eq.npp) k=npp-1
         dps=arg-psrad(k)
-        temp=spli(ct,npmx,k,dps)
+        temp=fspli(ct,npmx,k,dps)
       endif
       return
       end
@@ -3686,11 +3686,11 @@ c
       if(iprof.eq.0) then
         fzeff=zeff
       else
-        call locate(psrad,npp,ps,k)
+        call vlocate(psrad,npp,ps,k)
         if(k.eq.0) k=1
         if(k.eq.npp) k=npp-1
         dps=ps-psrad(k)
-        fzeff=spli(cz,npmx,k,dps)
+        fzeff=fspli(cz,npmx,k,dps)
       endif
       return
       end
@@ -3947,18 +3947,18 @@ c
             y0m=y0/1.0d2
             r0m=r0/1.0d2
             z0m=z0/1.0d2
-            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
+            write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
      .                    -1.0d0,zero,zero,zero,one
           end if
           if(j.eq.1.and.k.eq.1) then
-            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
+            write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
      .       zero,zero,zero,zero,
      .       zero,zero,zero,zero,one,zero,zero,
      .       zero,zero,one,zero,zero,zero,zero,one
             ddr110=dd
           end if  
           if(j.eq.nrayr.and.k.eq.1) then
-            write(17,99) zero,ddr110,dd,ddi
+            write(617,99) zero,ddr110,dd,ddi
           end if          
         end do
       end do
@@ -3967,7 +3967,7 @@ c
 c
       if(nrayr.gt.1) then
         iproj=0
-        nfilp=8
+        nfilp=608
         call projxyzt(iproj,nfilp)
       end if
 c
@@ -4106,12 +4106,12 @@ c
             y0m=y0/1.0d2
             r0m=r0/1.0d2
             z0m=z0/1.0d2
-            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
+            write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
      .                    -1.0d0,zero,zero,zero,one
           end if
            if(j.eq.1.and.k.eq.1) then
-            write(17,99) zero,zero,zero,zero
-            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
+            write(617,99) zero,zero,zero,zero
+            write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
      .       zero,zero,zero,zero,
      .       zero,zero,zero,zero,one,zero,zero,
      .       zero,zero,one,zero,zero,zero,zero,one
@@ -4123,7 +4123,7 @@ c
 c
       if(nrayr.gt.1) then
         iproj=0
-        nfilp=8
+        nfilp=608
         call projxyzt(iproj,nfilp)
       end if
 c
@@ -4206,12 +4206,12 @@ c
             y0m=y0/1.0d2
             r0m=r0/1.0d2
             z0m=z0/1.0d2
-            write(33,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
+            write(633,111) izero,j,k,zero,x0m,y0m,r0m,z0m,
      .                    -1.0d0,zero,zero,zero,one
           end if
           if(j.eq.1.and.k.eq.1) then
-            write(17,99) zero,zero,zero,zero
-            write(4,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
+            write(617,99) zero,zero,zero,zero
+            write(604,99) zero,r0m,z0m,atan2(y0m,x0m)*180.0d0/pi,
      .       zero,zero,zero,zero,
      .       zero,zero,zero,zero,one,zero,zero,
      .       zero,zero,one,zero,zero,zero,zero,one
@@ -4223,7 +4223,7 @@ c
 c
       if(nrayr.gt.1) then
         iproj=0
-        nfilp=8
+        nfilp=608
         call projxyzt(iproj,nfilp)
       end if
 c
@@ -4290,17 +4290,17 @@ c
 c
       dps=rpsi-rpstab(ip)
 c
-      areai=spli(carea,nintp,ip,dps)
-      voli=spli(cvol,nintp,ip,dps)
-      dervoli=splid(cvol,nintp,ip,dps)
-      rrii=spli(crri,nintp,ip,dps)
+      areai=fspli(carea,nintp,ip,dps)
+      voli=fspli(cvol,nintp,ip,dps)
+      dervoli=fsplid(cvol,nintp,ip,dps)
+      rrii=fspli(crri,nintp,ip,dps)
 c
       if(intp.eq.0) return
 c
-      rbavi=spli(crbav,nintp,ip,dps)
-      bmxi=spli(cbmx,nintp,ip,dps)
-      bmni=spli(cbmn,nintp,ip,dps)
-      fci=spli(cfc,nintp,ip,dps)
+      rbavi=fspli(crbav,nintp,ip,dps)
+      bmxi=fspli(cbmx,nintp,ip,dps)
+      bmni=fspli(cbmn,nintp,ip,dps)
+      fci=fspli(cfc,nintp,ip,dps)
 c
       return
       end
@@ -4319,9 +4319,9 @@ c
       if(ip.eq.0) ip=1
       if(ip.eq.nintp) ip=nintp-1
       dps=rpsi-rpstab(ip)
-      ratjai=spli(cratja,nintp,ip,dps)
-      ratjbi=spli(cratjb,nintp,ip,dps)
-      ratjpli=spli(cratjpl,nintp,ip,dps)
+      ratjai=fspli(cratja,nintp,ip,dps)
+      ratjbi=fspli(cratjb,nintp,ip,dps)
+      ratjpli=fspli(cratjpl,nintp,ip,dps)
       return
       end
 c
@@ -4488,7 +4488,7 @@ c      ratiovgr=2.0d0*anpr/derdnm*vgm
       alpha=2.0d0*akim*ratiovgr
       if(alpha.lt.0.0d0) then
         ierr=94
-        print*,'    IERR = ', ierr,'  alpha negative'
+        write(*,*) '    IERR = ', ierr,'  alpha negative'
       end if
 c
       ithn=1
@@ -4583,7 +4583,7 @@ c
 c
         if(dble(anpr2).lt.0.0d0.and.dimag(anpr2).lt.0.0d0) then
           anpr2=0.0d0
-          print*,'    Y =',yg,'   nperp2 < 0'
+          write(*,*) '    Y =',yg,'   nperp2 < 0'
 c          ierr=99
           go to 999
         end if
@@ -4593,7 +4593,7 @@ c
       end do
 c
   999 continue
-       if(i.gt.imx) print*,'    i>imx ',yg,errnpr,i
+       if(i.gt.imx) write(*,*) '    i>imx ',yg,errnpr,i
 c
       anpr=sqrt(anpr2)
       anprr=dble(anpr)
@@ -5513,7 +5513,7 @@ c     rzfc=(1+Zeff)/fc
         end do
         nhn=nhn-1
       CASE DEFAULT
-        print*,'ieccd undefined'
+        write(*,*) 'ieccd undefined'
        
       end select
 c
@@ -5913,7 +5913,7 @@ c     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
 c
       write(nfile,*) ' '
 c
-      write(12,99) istep,st,psinv11,rtimn,rtimx
+      write(612,99) istep,st,psinv11,rtimn,rtimx
       return
  99   format(i5,12(1x,e16.8e3))
 111   format(3i5,12(1x,e16.8e3))
@@ -6200,15 +6200,15 @@ c        call ratioj(rhopfi,ratjamx,ratjbmx,ratjplmx)
       pins_085=0.0d0
 
       xrhot=0.2d0
-      call locate(rhotv,nd,xrhot,i1)
+      call vlocate(rhotv,nd,xrhot,i1)
       call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
      .  xrhot,pins_02)
       xrhot=0.5d0
-      call locate(rhotv,nd,xrhot,i1)
+      call vlocate(rhotv,nd,xrhot,i1)
       call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
      .  xrhot,pins_05)
       xrhot=0.85d0
-      call locate(rhotv,nd,xrhot,i1)
+      call vlocate(rhotv,nd,xrhot,i1)
       call intlin(rhotv(i1),pins(i1),rhotv(i1+1),pins(i1+1),
      .  xrhot,pins_085)
      
@@ -6238,16 +6238,16 @@ c     dPdV [MW/m^3], Jcd [MA/m^2]
       if(ieccd.eq.0) currt=0.0d0
       currtka=currt*1.0d3
 
-      write(6,*)' '
-      write(6,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '//
+      write(*,*)' '
+      write(*,*)'#beta0 alpha0 Icd Pa Jphip dPdVp '//
      .'rhotj rhotjava rhotp rhotpav drhotjava drhotpav ratjamx '//
      .'ratjbmx stmx psipol chipol index_rt Jphimx dPdVmx drhotj drhotp'
-      write(6,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp,
+      write(*,99) beta0,alpha0,currtka,pabstot,ajphip,dpdvp,
      .  rhotjfi,rhotjava,rhotp,rhotpav,
      .  drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
      .  real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
      
-      write(7,99) currtka,pabstot,ajphip,dpdvp,
+      write(607,99) currtka,pabstot,ajphip,dpdvp,
      .  rhotjfi,rhotjava,rhotp,rhotpav,
      .  drhotjava,drhotpav,ratjamx,ratjbmx,stmx,psipol,chipol,
      .  real(index_rt),ajmxfi,dpdvmx,drhotjfi,drhotp
@@ -6262,7 +6262,7 @@ c     dPdV [MW/m^3], Jcd [MA/m^2]
         end if
         pinsr=0.0d0
         if(pabstot.gt.0) pinsr=pins(i)/pabstot
-        write(48,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
+        write(648,99) psin,rhotv(i),ajphiv(i),ajcdbv(i),dpdv(i),
      .   currins(i),pins(i),pinsr,real(index_rt)
       end do
 
@@ -6318,7 +6318,7 @@ c
         ypk=yy(2)
         rte1=0.0d0
         yye=ypk*emn1
-        call locate(yy,nd,yye,ie2)
+        call vlocate(yy,nd,yye,ie2)
         call intlin(yy(ie2),xx(ie2),yy(ie1+2),xx(ie1+2),yye,rte2)
       end if
 c
diff --git a/src/grayl.f b/src/grayl.f
index 0f0a9f6..3de8025 100644
--- a/src/grayl.f
+++ b/src/grayl.f
@@ -1,4 +1,4 @@
-      subroutine locate(xx,n,x,j)
+      subroutine vlocate(xx,n,x,j)
       implicit real*8(a-h,o-z)
       dimension xx(n)
 c
@@ -260,25 +260,25 @@ c
 c     spline routines: begin
 c
 c
-      function spli(cspli,n,k,dx)
+      function fspli(cspli,n,k,dx)
       implicit real*8(a-h,o-z)
       dimension cspli(n,4)
-      spli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4)))
+      fspli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4)))
       return
       end
 c
 c
 c
-      function splid(cspli,n,k,dx)
+      function fsplid(cspli,n,k,dx)
       implicit real*8(a-h,o-z)
       dimension cspli(n,4)
-      splid=cspli(k,2)+dx*(2.0d0*cspli(k,3)+3.0d0*dx*cspli(k,4))
+      fsplid=cspli(k,2)+dx*(2.0d0*cspli(k,3)+3.0d0*dx*cspli(k,4))
       return
       end
 c
 c
 c
-      subroutine difcs(x,y,n,iopt,c,ier)
+      subroutine difcsg(x,y,n,iopt,c,ier)
       implicit real*8(a-h,o-z)
       dimension x(n),y(n),c(n*4)
       jmp  =1
@@ -383,7 +383,7 @@ c
 c
       subroutine difcsn(xx,yy,nmx,n,iopt,cc,ier)
 c
-c     same as difcs but with dimension(xx,yy) = nmx > n
+c     same as difcsg but with dimension(xx,yy) = nmx > n
 c
       implicit real*8(a-h,o-z)
 c
@@ -689,7 +689,7 @@ c
       t(1)=acos(x)
       h(1)=tau*t(1)
       b0=besi0(h(1))
-      b1=besi1(h(1))
+      b1=besi1g(h(1))
       z=-1.0d0
       go to 3
 c
@@ -1220,7 +1220,7 @@ c
       e=.true.
       go to 2
 c
-      entry besi1(x)
+      entry besi1g(x)
 c
       e=.false.
     2 l=.true.
@@ -8485,13 +8485,13 @@ c  error codes and messages.
 c
 c
 c
-      subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
-c  subroutine splder evaluates in a number of points x(i),i=1,2,...,m
+      subroutine gsplder(t,n,c,k,nu,x,y,m,wrk,ier)
+c  subroutine gsplder evaluates in a number of points x(i),i=1,2,...,m
 c  the derivative of order nu of a spline s(x) of degree k,given in
 c  its b-spline representation.
 c
 c  calling sequence:
-c     call splder(t,n,c,k,nu,x,y,m,wrk,ier)
+c     call gsplder(t,n,c,k,nu,x,y,m,wrk,ier)
 c
 c  input parameters:
 c    t    : array,length n, which contains the position of the knots.
@@ -8907,15 +8907,15 @@ c  the zeros of s(x) are arranged in increasing order.
       end
 
 c
-      subroutine profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
-c  if iopt=0 subroutine profil calculates the b-spline coefficients of
+      subroutine sprofil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
+c  if iopt=0 subroutine sprofil calculates the b-spline coefficients of
 c  the univariate spline f(y) = s(u,y) with s(x,y) a bivariate spline of
 c  degrees kx and ky, given in the b-spline representation.
 c  if iopt = 1 it calculates the b-spline coefficients of the univariate
 c  spline g(x) = s(x,u)
 c
 c  calling sequence:
-c     call profil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
+c     call sprofil(iopt,tx,nx,ty,ny,c,kx,ky,u,nu,cu,ier)
 c
 c  input parameters:
 c   iopt  : integer flag, specifying whether the profile f(y) (iopt=0)