From 2c90c5f2cfe688890a7dcf75d830657ea2d4fd98 Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Thu, 19 Jun 2014 08:13:52 +0000
Subject: [PATCH] renamed clashing common block names; fixed wrong calling of
 splining routines; fixed spline evaluation functions to work with uneven x
 grid spacing; removed debugging prints

---
 Makefile             |   2 +-
 src/gray-externals.f | 196 ++++++++++++++++++++++---------------------
 src/gray_main.f90    |  12 ---
 src/grayl.f          |   2 +-
 4 files changed, 102 insertions(+), 110 deletions(-)

diff --git a/Makefile b/Makefile
index 97f3803..3b5a4c4 100644
--- a/Makefile
+++ b/Makefile
@@ -97,7 +97,7 @@ clean:
 # Dependencies
 # ------------
 gray_main.o: const_and_precisions.o
-gray-externals.o: green_func_p.o reflections.o beamdata.o
+gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
 green_func_p.o: const_and_precisions.o
 scatterspl.o: const_and_precisions.o
 beamdata.o: const_and_precisions.o
diff --git a/src/gray-externals.f b/src/gray-externals.f
index 5d679f1..4b03a35 100644
--- a/src/gray-externals.f
+++ b/src/gray-externals.f
@@ -114,7 +114,7 @@ c
       common/ist/istpr0,istpl0
       common/istgr/istpr,istpl
 c
-      common/psival/psinv
+      common/psinv/psinv
       common/psinv11/psinv11
       common/ierr/ierr
       common/taumnx/taumn,taumx,pabstot,currtot
@@ -134,7 +134,7 @@ c
       common/index_rt/index_rt
       common/ipass/ipass
       common/rwallm/rwallm
-      common/bound/zbmin,zbmax
+      common/zbound/zbmin,zbmax
 
       pabstot=0.0d0
       currtot=0.0d0
@@ -306,7 +306,7 @@ c
       common/btot/btot
       common/xgxg/xg
       common/ygyg/yg
-      common/dens/dens,ddens
+      common/dddens/dens,ddens
       common/tete/tekev
       common/absor/alpha,effjcd,akim,tau0
       common/densbnd/psdbnd
@@ -732,9 +732,6 @@ c      anz0c=-cos(cvdr*beta0)*sin(cvdr*alpha0)
  
       anx0c=(anr0c*x00-anphi0c*y00)/r00
       any0c=(anr0c*y00+anphi0c*x00)/r00
-c
-      print*,'    input file read'
-!      call myflush
 c
 c     read data for Te , ne , Zeff from file if iprof>0
 c
@@ -748,8 +745,6 @@ c     read profiles from input arguments
         call profdata(nrho, psijet, te, dne, zeff)
 c        close(nprof)
       end if
-      print*,'    profiles fitted'
-!     call myflush
 c
 c     read equilibrium data  from file if iequil=2
 c
@@ -762,8 +757,6 @@ c     .        status= 'unknown', unit=neqdsk)
         call equidata(ijetto, mr, mz, r, z, psin, psiax, psibnd, 
      .      rax, zax, nbnd, rbnd, zbnd, nrho, psijet, f, qsf)
 c        close(neqdsk)
-      print*,'    equilibrium fitted'
-!     call myflush
         
 c     print density, temperature, safecty factor, toroidal current dens
 c     versus psi, rhop, rhot
@@ -964,16 +957,16 @@ c
       if(alpha0.gt.alphastv(1).and.alpha0.lt.alphastv(nisteer)) then
         call vlocate(alphastv,nisteer,alpha0,k)
         dal=alpha0-alphastv(k)
-        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)
+        betst=fspli(cbeta,nstrmx,k,dal)
+        x00=fspli(cx0,nstrmx,k,dal)
+        y00=fspli(cy0,nstrmx,k,dal)
+        z00=fspli(cz0,nstrmx,k,dal)
+        wcsi=fspli(cwaist1,nstrmx,k,dal)
+        weta=fspli(cwaist2,nstrmx,k,dal)
+        rcicsi=fspli(crci1,nstrmx,k,dal)
+        rcieta=fspli(crci2,nstrmx,k,dal)
+        phiw=fspli(cphi1,nstrmx,k,dal)
+        phir=fspli(cphi2,nstrmx,k,dal)
       else
         write(*,*) '    alpha0 outside table range !!!'
         if(alpha0.ge.alphastv(nisteer)) ii=nisteer
@@ -1043,7 +1036,7 @@ c
       common/ipsn/ipsinorm
       common/sspl/sspl
       common/nfile/neqdsk,nprof
-      common/bound/zbmin,zbmax
+      common/zbound/zbmin,zbmax
       common/sgnib/sgnbphi,sgniphi
       common/factb/factb
       common/ixp/ixp
@@ -1106,7 +1099,7 @@ c     psi function
       psia0=psiedge-psiaxis
       psia=psia0*factb
       sgniphi=sign(1.0d0,-psia0)
-cc
+c
 c      do j=1,nz
 c        do i=1,nr
 c          write(620,2021) rv(i),zv(j),psin(i,j)
@@ -1172,21 +1165,21 @@ c
         nsr=nr/4+4
         nsz=nz/4+4
         call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
-     .      rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier)
+     .      rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier)
 c       if ier=-1 data are fitted using sspl=0    
         if(ier.eq.-1) then
           sspl=0.0d0
           nsr=nr/4+4
           nsz=nz/4+4
           call scatterspl(rv1d,zv1d,fvpsi,wpsi,nrz,kspl,sspl,
-     .        rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,coeff,ier)
+     .        rmnm,rmxm,zmnm,zmxm,tr,nsr,tz,nsz,cc,ier)
         end if
         nsrt=nsr
         nszt=nsz
       end if
-cc
-cc     re-evaluate psi on the grid using the spline (only for debug and cniteq)
-cc
+c
+c     re-evaluate psi on the grid using the spline (only for debug and cniteq)
+c
 c      call dierckx_bispev(tr,nsr,tz,nsz,cc,kspl,kspl,rv,nr,zv,nz,ffvpsi,
 c     .   wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ier)
 c
@@ -1354,15 +1347,15 @@ c
 c     compute B_toroidal on axis
 
       btaxis=fpol(1)/rmaxis
-      btrcen=abs(btrcen)*factb
+      btrcen=btrcen*factb
       write(*,'(a,f8.4)') 'factb =   ',factb
-      write(*,'(a,f8.4)') '|BT_centr|= ',btrcen
+      write(*,'(a,f8.4)') '|BT_centr|= ',abs(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)
+      rrtor=sqrt(abs(phitedge/btrcen)/pi)
       call rhopol
 c      write(*,*) rhotsx,phitedge,rrtor,abs(psia)
 
@@ -1417,7 +1410,7 @@ c
       parameter(n=2,ldfjac=n,lwa=(n*(n+13))/2)
       dimension xvec(n),fvec(n),fjac(ldfjac,n),wa(lwa)
       external fcnox
-      common/psival/psinv
+      common/psinv/psinv
       xvec(1)=rz
       xvec(2)=zz
       tol = sqrt(comp_eps)
@@ -1479,7 +1472,7 @@ c
       implicit real*8 (a-h,o-z)
       dimension x(n),fvec(n),fjac(ldfjac,n)
       common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
-      common/psival/psinv
+      common/psinv/psinv
       common/cnpsi/h
       common/pareq1/psia
       call equinum(x(1),x(2))
@@ -1505,7 +1498,7 @@ c
       common/psinr/psinr
       common/qpsi/qpsi
       common/eqnn/nr,nz,nrho,npp,nintp
-      common/dens/dens,ddens
+      common/dddens/dens,ddens
 c
       write(645,*) ' #psi rhot ne Te q Jphi'
       psin=0.0d0
@@ -1739,14 +1732,14 @@ c
         rhotnr(k+1)=rhotnr(k)+drhot
       end do
       rhotsx=rhotnr(nr)
-      do k=1,nr
+      do k=2,nr
         rhotnr(k)=sqrt(rhotnr(k)/rhotnr(nr))
       end do
 c
 c     spline interpolation of rhotor
 c
       iopt=0
-      call difcsg(psinr,rhotnr,nr,iopt,crhot,ier)
+      call difcsn(psinr,rhotnr,nnw,nr,iopt,crhot,ier)
       return
       end
 
@@ -1757,11 +1750,10 @@ c
       common/psinr/psinr
       common/eqnn/nr,nz,nrho,npp,nintp
       common/cq/cq
-      irt=int((nrho-1)*psi+1)
-      if(irt.eq.0) irt=1
-      if(irt.eq.nrho) irt=nrho-1
+      call vlocate(psinr,nrho,psi,irt)
+      irt=max(1,min(irt,nrho-1))
       dps=psi-psinr(irt)
-      fq_eq=fspli(cq,nrho,irt,dps)
+      fq_eq=fspli(cq,nnw,irt,dps)
       return
       end
 
@@ -1773,11 +1765,10 @@ c
       common/eqnn/nr,nz,nrho,npp,nintp
       common/crhot/crhot
 c
-      irt=int((nrho-1)*psi+1)
-      if(irt.eq.0) irt=1
-      if(irt.eq.nrho) irt=nrho-1
+      call vlocate(psinr,nrho,psi,irt)
+      irt=max(1,min(irt,nrho-1))
       dps=psi-psinr(irt)
-      frhotor_eq=fspli(crhot,nrho,irt,dps)
+      frhotor_eq=fspli(crhot,nnw,irt,dps)
       return
       end
 
@@ -1796,11 +1787,13 @@ c
       common/crhotq/crhotq
 
       rpsi=sqrt(psi)
-      ip=int((nintp-1)*rpsi+1)
-      if(ip.eq.0) ip=1
-      if(ip.eq.nintp) ip=nintp-1
+c      ip=int((nintp-1)*rpsi+1)
+c      if(ip.eq.0) ip=1
+c      if(ip.eq.nintp) ip=nintp-1
+      call vlocate(rpstab,nintp,rpsi,ip)
+      ip=max(1,min(ip,nintp-1))
       dps=rpsi-rpstab(ip)
-      frhotor_av=fspli(crhotq,nintp,ip,dps)
+      frhotor_av=fspli(crhotq,nnintp,ip,dps)
       return
       end
       
@@ -1816,14 +1809,18 @@ c
       common/coffrp/crp
      
       dr=1.0d0/dble(nnr-1)
-      do i=1,nnr
+      do i=2,nnr-1
         rhop(i)=(i-1)*dr
         psin=rhop(i)*rhop(i)
-        rhot(i)=frhotor(psin)
+        rhot(i)=frhotor_eq(psin)
         wp(i)=1.0d0
       end do 
-        wp(1)=1.0d3
-        wp(nnr)=1.0d3
+      rhop(1)=0.0d0
+      rhot(1)=0.0d0
+      wp(1)=1.0d3
+      rhop(nnr)=1.0d0
+      rhot(nnr)=1.0d0
+      wp(nnr)=1.0d3
       
 c     spline interpolation of rhopol versus rhotor
       iopt=0
@@ -2156,7 +2153,7 @@ c
       common/cratj/cratja,cratjb,cratjpl
       common/crhotq/crhotq
       common/cnt/rup,zup,rlw,zlw
-      common/bound/zbmin,zbmax
+      common/zbound/zbmin,zbmax
       common/rarea/rarea
       common/phitedge/phitedge
       common/cdadrhot/cdadrhot
@@ -2478,11 +2475,13 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       common/pstab/rpstab
       common/eqnn/nr,nz,nrho,npp,nintp
       common/cdadrhot/cdadrhot
-      ip=int((nintp-1)*rpsi+1)
-      if(ip.eq.0) ip=1
-      if(ip.eq.nintp) ip=nintp-1
+c      ip=int((nintp-1)*rpsi+1)
+c      if(ip.eq.0) ip=1
+c      if(ip.eq.nintp) ip=nintp-1
+      call vlocate(rpstab,nintp,rpsi,ip)
+      ip=max(1,min(ip,nintp-1))
       dps=rpsi-rpstab(ip)
-      fdadrhot=fspli(cdadrhot,nintp,ip,dps)
+      fdadrhot=fspli(cdadrhot,nnintp,ip,dps)
       return
       end
 
@@ -2493,11 +2492,13 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       common/pstab/rpstab
       common/eqnn/nr,nz,nrho,npp,nintp
       common/cdvdrhot/cdvdrhot
-      ip=int((nintp-1)*rpsi+1)
-      if(ip.eq.0) ip=1
-      if(ip.eq.nintp) ip=nintp-1
+c      ip=int((nintp-1)*rpsi+1)
+c      if(ip.eq.0) ip=1
+c      if(ip.eq.nintp) ip=nintp-1
+      call vlocate(rpstab,nintp,rpsi,ip)
+      ip=max(1,min(ip,nintp-1))
       dps=rpsi-rpstab(ip)
-      fdvdrhot=fspli(cdvdrhot,nintp,ip,dps)
+      fdvdrhot=fspli(cdvdrhot,nnintp,ip,dps)
       return
       end
 
@@ -3140,7 +3141,7 @@ c
       common/dxgyg/derxg,deryg
       common/iieq/iequil
       common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
-      common/psival/psinv
+      common/psinv/psinv
       common/sgnib/sgnbphi,sgniphi
 c
       xg=0.0d0
@@ -3292,14 +3293,14 @@ c
 c
       common/parqq/q0,qa,alq
       common/parban/b0,rr0m,zr0m,rpam
-      common/psival/psinv
+      common/psinv/psinv
       common/pareq1/psia
       common/densbnd/psdbnd
       common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
       common/xgxg/xg
       common/dxgdps/dxgdpsi
       common/xgcn/xgcn
-      common/dens/dens,ddens
+      common/dddens/dens,ddens
       common/sgnib/sgnbphi,sgniphi
       common/bmxmn/bmxi,bmni
       common/fc/fci
@@ -3381,7 +3382,7 @@ c
       dimension tfp(nrest),cfp(nrest),wrkfd(nrest)
 c
       common/eqnn/nr,nz,nrho,npp,nintp
-      common/psival/psinv
+      common/psinv/psinv
       common/derip/dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,ffpv
       common/pareq1/psia
       common/pareq2/btaxis,rmaxis,zmaxis,rmnm,rmxm,zmnm,zmxm,dr,dz
@@ -3415,7 +3416,7 @@ c
       zzs(1)=zpsim
       nsr=nsrt
       nsz=nszt
-      call fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
+      call dierckx_fpbisp(tr,nsr,tz,nsz,ccspl,3,3,
      .     rrs,nrs,zzs,nzs,ffspl,wrk(1),wrk(5),iwrk(1),iwrk(2))
       psinv=(ffspl(1)-psinop)/psiant
       if(psinv.lt.0.0d0)
@@ -3427,7 +3428,8 @@ c
       kkz=3-nuz
       iwr=1+(nr-nur-4)*(nz-nuz-4)
       iwz=iwr+4-nur
-      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,kkr,kkz,
+      call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc10,
+     .     kkr,kkz,
      .     rrs,nrs,zzs,nzs,ffspl,cc10(iwr),cc10(iwz),iwrk(1),iwrk(2))
       dpsidr= ffspl(1)*psia
 c
@@ -3437,7 +3439,8 @@ c
       kkz=3-nuz
       iwr=1+(nr-nur-4)*(nz-nuz-4)
       iwz=iwr+4-nur
-      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,kkr,kkz,
+      call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc01,
+     .     kkr,kkz,
      .     rrs,nrs,zzs,nzs,ffspl,cc01(iwr),cc01(iwz),iwrk(1),iwrk(2))
       dpsidz= ffspl(1)*psia
 c
@@ -3447,7 +3450,8 @@ c
       kkz=3-nuz
       iwr=1+(nr-nur-4)*(nz-nuz-4)
       iwz=iwr+4-nur
-      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,kkr,kkz,
+      call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc20,
+     .     kkr,kkz,
      .     rrs,nrs,zzs,nzs,ffspl,cc20(iwr),cc20(iwz),iwrk(1),iwrk(2))
       ddpsidrr= ffspl(1)*psia
 c
@@ -3457,7 +3461,8 @@ c
       kkz=3-nuz
       iwr=1+(nr-nur-4)*(nz-nuz-4)
       iwz=iwr+4-nur
-      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,kkr,kkz,
+      call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc02,
+     .     kkr,kkz,
      .     rrs,nrs,zzs,nzs,ffspl,cc02(iwr),cc02(iwz),iwrk(1),iwrk(2))
       ddpsidzz= ffspl(1)*psia
 c
@@ -3467,7 +3472,8 @@ c
       kkz=3-nuz
       iwr=1+(nr-nur-4)*(nz-nuz-4)
       iwz=iwr+4-nur
-      call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,kkr,kkz,
+      call dierckx_fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc11,
+     .     kkr,kkz,
      .     rrs,nrs,zzs,nzs,ffspl,cc11(iwr),cc11(iwz),iwrk(1),iwrk(2))
       ddpsidrz= ffspl(1)*psia
 c
@@ -3548,12 +3554,12 @@ c
 c
       subroutine sub_xg_derxg
       implicit real*8 (a-h,o-z)
-      common/psival/psinv
+      common/psinv/psinv
       common/pareq1/psia
       common/xgxg/xg
       common/dxgdps/dxgdpsi
       common/xgcn/xgcn
-      common/dens/dens,ddenspsin
+      common/dddens/dens,ddenspsin
       xg=0.0d0
       dxgdpsi=0.0d0
 c      if(psinv.le.psdbnd.and.psinv.ge.0) then
@@ -3576,7 +3582,7 @@ c
       common/denspp/psnpp,aad,bbd,ccd
       common/iipr/iprof
       common/pardens/dens0,aln1,aln2
-      common/dens/dens,ddens
+      common/dddens/dens,ddens
       common/coffdt/tfn
       common/coffdnst/nsfd
       common/cofffn/cfn
@@ -3644,8 +3650,7 @@ c
         temperature=(te0-dte0)*proft+dte0
       else
         call vlocate(psrad,npp,arg,k)
-        if(k.eq.0) k=1
-        if(k.eq.npp) k=npp-1
+        k=max(1,min(k,npp-1))
         dps=arg-psrad(k)
         temperature=fspli(ct,npmx,k,dps)
       endif
@@ -3672,8 +3677,7 @@ c
         fzeff=zeff
       else
         call vlocate(psrad,npp,ps,k)
-        if(k.eq.0) k=1
-        if(k.eq.npp) k=npp-1
+        k=max(1,min(k,npp-1))
         dps=ps-psrad(k)
         fzeff=fspli(cz,npmx,k,dps)
       endif
@@ -4238,17 +4242,17 @@ c
 c
       dps=rpsi-rpstab(ip)
 c
-      areai=fspli(carea,nintp,ip,dps)
-      voli=fspli(cvol,nintp,ip,dps)
-      dervoli=fsplid(cvol,nintp,ip,dps)
-      rrii=fspli(crri,nintp,ip,dps)
+      areai=fspli(carea,nnintp,ip,dps)
+      voli=fspli(cvol,nnintp,ip,dps)
+      dervoli=fsplid(cvol,nnintp,ip,dps)
+      rrii=fspli(crri,nnintp,ip,dps)
 c
       if(intp.eq.0) return
 c
-      rbavi=fspli(crbav,nintp,ip,dps)
-      bmxi=fspli(cbmx,nintp,ip,dps)
-      bmni=fspli(cbmn,nintp,ip,dps)
-      fci=fspli(cfc,nintp,ip,dps)
+      rbavi=fspli(crbav,nnintp,ip,dps)
+      bmxi=fspli(cbmx,nnintp,ip,dps)
+      bmni=fspli(cbmn,nnintp,ip,dps)
+      fci=fspli(cfc,nnintp,ip,dps)
 c
       return
       end
@@ -4267,9 +4271,9 @@ c
       if(ip.eq.0) ip=1
       if(ip.eq.nintp) ip=nintp-1
       dps=rpsi-rpstab(ip)
-      ratjai=fspli(cratja,nintp,ip,dps)
-      ratjbi=fspli(cratjb,nintp,ip,dps)
-      ratjpli=fspli(cratjpl,nintp,ip,dps)
+      ratjai=fspli(cratja,nnintp,ip,dps)
+      ratjbi=fspli(cratjb,nnintp,ip,dps)
+      ratjpli=fspli(cratjpl,nnintp,ip,dps)
       return
       end
 c
@@ -4290,7 +4294,7 @@ c
       common/parban/b0,rr0m,zr0m,rpam
       common/absor/alpha,effjcd,akim,tau0
 c
-      common/psival/psinv
+      common/psinv/psinv
       common/sgnib/sgnbphi,sgniphi
       common/bmxmn/bmxi,bmni
       common/fc/fci
@@ -4364,7 +4368,7 @@ c
 c
       common/absor/alpha,effjcd,akim,tau
 c
-      common/psival/psinv
+      common/psinv/psinv
       common/tete/tekev
       common/amut/amu
       common/zz/Zeff
@@ -5381,7 +5385,7 @@ c
       
       common/ieccd/ieccd
       common/tete/tekev
-      common/dens/dens,ddens
+      common/dddens/dens,ddens
       common/zz/Zeff
       common/btot/btot
       common/bmxmn/bmax,bmin
@@ -5713,7 +5717,7 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       common/nplr/anpl,anpr
       common/fc/fc
       common/ncl/hb
-      common/psival/psinv
+      common/psinv/psinv
       common/amut/amu
       common/tete/tekev
       common/zz/Zeff
@@ -5747,7 +5751,7 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       parameter(kwrk=nnintp+nlam+njest+nlest+3)
       parameter(lw01=nnintp*4+nlam*3+nnintp*nlam)
       
-      external fpbisp
+      external dierckx_fpbisp
 
       dimension xxs(1),yys(1),ffs(1)
       dimension ch01(lw01),ch((njest-4)*(nlest-4))
@@ -5764,13 +5768,13 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       xxs(1)=sqrt(psi)
       yys(1)=alam
       
-      call fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs,
+      call dierckx_fpbisp(tjp,njp,tlm,nlm,ch,ksp,ksp,xxs,1,yys,1,ffs,
      .            wrk(1),wrk(5),iwrk(1),iwrk(2))
       fv=ffs(1)
 
       iwp=1+(njp-4)*(nlm-5)
       iwl=iwp+4
-      call fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2,
+      call dierckx_fpbisp(tjp(1),njp,tlm(2),nlm-2,ch01,3,2,
      .     xxs,1,yys,1,ffs,ch01(iwp),ch01(iwl),iwrk(1),iwrk(2))
       dfv=ffs(1)
 
@@ -6580,7 +6584,7 @@ c     wave vector and electric field after reflection in lab frame
       common/limiter/rlim,zlim,nlim
       common/anv/anv
       common/dsds/dst      
-      common/psival/psinv
+      common/psinv/psinv
       common/densbnd/psdbnd
       common/dstvac/dstvac
 c     ivac=1 first interface plasma-vacuum
diff --git a/src/gray_main.f90 b/src/gray_main.f90
index 1b31b4d..9165fcf 100644
--- a/src/gray_main.f90
+++ b/src/gray_main.f90
@@ -33,21 +33,11 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
       
 ! read data plus initialization
   index_rt=1
-  print*,'GRAY started'
-!  call myflush
   call prfile
-  print*,'  file headers written'
-!  call myflush
   call paraminit
-  print*,'  variables initialized'
-!  call myflush
   call read_data(ijetto, mr, mz, r, z, psin, psiax, psibnd, rax, zax, &
     nbnd, rbnd, zbnd, nrho, psijet, f, te, dne, zeff, qsf, powin)
-  print*,'  spline computed'
-!  call myflush
   call vectinit
-  print*,'  beam arrays allocated'
-!  call myflush
   if(iercom.eq.0) then
     if(igrad.eq.0) call ic_rt
     if(igrad.gt.0) call ic_gb
@@ -57,8 +47,6 @@ subroutine gray_main(ijetto, mr, mz, r, z, psin, psiax, psibnd, &
     write(*,*) '    IERR = ', ierr
     return
   end if
-  print*,'  initial conditions set'
-!  call myflush
 
 ! beam/ray propagation
   call gray_integration
diff --git a/src/grayl.f b/src/grayl.f
index 4ea42f5..55762d9 100644
--- a/src/grayl.f
+++ b/src/grayl.f
@@ -7733,7 +7733,7 @@ c  we partition the working space and evaluate the partial derivative
       end
 
 
-      subroutine dierckx_coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,
+      subroutine coeff_parder(tx,nx,ty,ny,c,kx,ky,nux,nuy,
      * wrk,lwrk,ier)
 c  subroutine parder evaluates on a grid (x(i),y(j)),i=1,...,mx; j=1,...
 c  ,my the partial derivative ( order nux,nuy) of a bivariate spline