diff --git a/src/gray.f b/src/gray.f
index fd4a4f6..ae70b6c 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -5686,7 +5686,11 @@ c local constants
 c local variables
       integer :: lrm,ithn
       real(wp_) :: amu,ratiovgr,rbn,rbx,zeff
+      integer :: npar
+      real(wp_) :: cst2
+      real(wp_), dimension(:), allocatable :: eccdpar
 c common/external functions/variables
+      real(wp_), external :: fjch,fjncl,fjch0
       integer :: nhmin,nhmax,iokhawa,ierr
       real(wp_) :: xg,yg,anpl,anpr,vgm,derdnm,ak0,akinv,fhz,sox,
      . anprre,anprim,alpha,effjcd,akim,tau,psinv,tekev,dens,ddens,
@@ -5712,6 +5716,16 @@ c
       common/btot/btot
       common/bmxmn/bmax,bmin
       common/fc/fc
+c
+      interface
+        subroutine setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,
+     .    cst2,eccdpar,npar)
+        use const_and_precisions, only : wp_
+        integer :: ieccd,npar
+        real(wp_) :: amu,Zeff,rbn,rbx,fc,psinv,cst2
+        real(wp_), dimension(:), allocatable :: eccdpar
+        end subroutine setcdcoeff
+      end interface
 c
 c     absorption computation
 c
@@ -5746,52 +5760,56 @@ c
         zeff=fzeff(psinv)
         rbn=btot/bmin
         rbx=btot/bmax
-        call eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez,
-     *  dens,psinv,ieccd,nhmin,nhmax,ithn,effjcd,iokhawa,ierr)
-      end if
 
+        call setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,cst2,
+     .      eccdpar,npar)
+        select case(ieccd) 
+        case(1)
+c         cohen model      
+          call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
+     .       ithn,cst2,fjch,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
+        case(2:9)
+c         no trapping
+          call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
+     .       ithn,cst2,fjch0,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
+        case(10:11)
+c         neoclassical model
+          call eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmin,nhmax,
+     .       ithn,cst2,fjncl,eccdpar(1:npar),npar,effjcd,iokhawa,ierr)
+        CASE DEFAULT
+          effjcd=0.0_wp_
+          print*,'ieccd undefined'
+          ierr=89
+          return
+        end select
+c
+        deallocate(eccdpar)
+      end if
+c
       return
       end
 c
 c
 c
-      subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez,
-     .  dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr)
-      use const_and_precisions, only : wp_,pi,qesi=>e_,mesi=>me_,
-     . vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
+    
+      subroutine setcdcoeff(amu,zeff,rbn,rbx,fc,psinv,ieccd,
+     .    cst2,eccdpar,npar)
+      use const_and_precisions, only : wp_
       use green_func_p, only: SpitzFuncCoeff
       use conical, only : fconic
       use magsurf_data, only : ch,tjp,tlm,njpt,nlmt
       use dierckx, only : profil
-      implicit none
 c local constants
-      real(wp_), parameter :: mc2m2=1.0_wp_/mc2**2,
-     . canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi)
+      real(wp_), parameter :: cst2min=1.0e-6_wp_
       integer, parameter :: ksp=3
 c arguments
-      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
+      integer :: ieccd,npar
+      real(wp_) :: amu,Zeff,rbn,rbx,fc,psinv,cst2
+      real(wp_), dimension(:), allocatable :: eccdpar ! dimension(max(5,3+nlmt))
 c local variables
-      integer :: nhn,njp,nlm,npar
-      real(wp_) :: anum,denom,resp,resj,coullog,anucc,alams,
-     . fp0s,pa,cst2
+      integer :: njp,nlm
+      real(wp_) :: alams,fp0s,pa
       real(wp_) :: chlm(nlmt)
-      real(wp_), dimension(3+2*nlmt) :: eccdpar ! dimension(max(5,3+nlmt))
-c common/external functions/variables
-      real(wp_), external :: fjch,fjncl,fjch0
-c
-      anum=0.0_wp_
-      denom=0.0_wp_
-      effjcd=0.0_wp_
-c
-      coullog=48.0_wp_-log(1.0e7_wp_*dens*mc2m2*amu**2)
-      anucc=canucc*dens*coullog
-c
-c      nhmx=nhmn
-c
-      eccdpar(1)=zeff
 c
       select case(ieccd) 
 c
@@ -5802,114 +5820,115 @@ c     rbx=B/B_max
 c     cst2=1.0_wp_-B/B_max
 c     alams=sqrt(1-B_min/B_max)
 c     Zeff < 31 !!!
-c     fp0s= P_a (alams)   
+c     fp0s= P_a (alams)
         cst2=1.0_wp_-rbx
-        if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_
+        if(cst2.lt.cst2min) cst2=0.0_wp_
         alams=sqrt(1.0_wp_-rbx/rbn)
         pa=sqrt(32.0_wp_/(Zeff+1.0_wp_)-1.0_wp_)/2.0_wp_
         fp0s=fconic(alams,pa,0)
+        npar=5
+        allocate(eccdpar(npar))
+        eccdpar(1)=zeff
         eccdpar(2)=rbn
         eccdpar(3)=alams
         eccdpar(4)=pa
         eccdpar(5)=fp0s
-        do nhn=nhmn,nhmx
-          call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
-     .                  fjch,eccdpar,5,resj,resp,iokhawa,ierr)
-          anum=anum+resj
-          denom=denom+resp
-        end do
-       
+c
       case(2:9)
         cst2=0.0_wp_
-        do nhn=nhmn,nhmx
-          call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
-     .                  fjch0,eccdpar,1,resj,resp,iokhawa,ierr)
-          anum=anum+resj
-          denom=denom+resp
-        end do
-     
+        npar=1
+        allocate(eccdpar(npar))
+        eccdpar(1)=zeff
+c
       case(10:11)
 c     neoclassical model:
 c     ft=1-fc trapped particle fraction
 c     rzfc=(1+Zeff)/fc     
         cst2=1.0_wp_-rbx
-        if(cst2.lt.1.0e-6_wp_) cst2=0.0_wp_
+        if(cst2.lt.cst2min) cst2=0.0_wp_
         call SpitzFuncCoeff(amu,Zeff,fc)
-        eccdpar(2) = fc
-        eccdpar(3) = rbx
-
         njp=njpt
         nlm=nlmt
         call profil(0,tjp,njp,tlm,nlm,ch,ksp,ksp,sqrt(psinv),nlm,chlm,
      .              ierr)
         if(ierr.gt.0) print*,'    Hlambda profil =',ierr
         npar=3+2*nlm
+        allocate(eccdpar(npar))
+        eccdpar(1)=zeff
+        eccdpar(2) = fc
+        eccdpar(3) = rbx
         eccdpar(4:3+nlm) = tlm
         eccdpar(4+nlm:npar) = chlm
-        do nhn=nhmn,nhmx
-          call curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
-     .                  fjncl,eccdpar,npar,resj,resp,iokhawa,ierr)
-          anum=anum+resj
-          denom=denom+resp
-        end do
-       
-       CASE DEFAULT
-       print*,'ieccd undefined'
-       
+c
       end select
 c
-c     effjpl = <J_parallel>/<p_d> /(B_min/<B>)  [A m /W]
-c
-      if(denom.gt.0.0_wp_) effjcd=-ceff*anum/(anucc*denom)
       return
       end
 c
 c
 c
-      subroutine curr_int(yg,anpl,anprre,amu,ex,ey,ez,nhn,ithn,cst2,
-     *   fcur,eccdpar,necp,resj,resp,iokhawa,ierr)
-      use const_and_precisions, only : wp_
+      subroutine eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmn,nhmx,
+     *   ithn,cst2,fcur,eccdpar,necp,effjcd,iokhawa,ierr)
+      use const_and_precisions, only : wp_,pi,qesi=>e_,mesi=>me_,
+     . vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
       use quadpack, only : dqagsmv
       implicit none
 c local constants
+      real(wp_), parameter :: mc2m2=1.0_wp_/mc2**2,
+     . canucc=2.0e13_wp_*pi*qe**4/(me**2*vc**3),ceff=qesi/(mesi*vcsi)
       real(wp_), parameter :: epsa=0.0_wp_,epsr=1.0e-2_wp_,xxcr=16.0_wp_
-      integer, parameter :: lw=5000,liw=lw/4,nfpp=15
+      real(wp_), parameter :: dumin=1.0e-6_wp_
+      integer, parameter :: lw=5000,liw=lw/4,nfpp=13
 c arguments
-      integer :: i,nhn,ithn,necp,iokhawa,ierr
-      real(wp_) :: yg,anpl,anprre,amu,cst2,resj,resp
+      integer :: i,nhmn,nhmx,ithn,necp,iokhawa,ierr
+      real(wp_) :: yg,anpl,anprre,dens,amu,cst2,effjcd
       real(wp_), dimension(necp) :: eccdpar
       complex(wp_) :: ex,ey,ez
 c local variables
-      integer :: neval,ier,last,npar
+      integer :: nhn,neval,ier,last,npar
       integer, dimension(liw) :: iw
-      real(wp_) :: anpl2,dnl,ygn,ygn2,resj1,resj2,rdu2,upltp,upltm,
-     . rdu,rdut,rdu2t,duu,uu1,uu2,xx1,xx2,ej,ej1,ej2,epp,uplp,uplm,
-     . cst
+      real(wp_) :: anpl2,dnl,ygn,ygn2,resji,rdu2,upltp,upltm,uplp,uplm,
+     . rdu,rdu2t,duu,uu1,uu2,xx1,xx2,resj,resp,eji,epp,anum,denom,
+     . cstrdut,anucc
       real(wp_), dimension(lw) :: w
       real(wp_), dimension(nfpp+necp) :: apar
+      real(wp_), dimension(0:1) :: uleft,uright
 c common/external functions/variables
       real(wp_), external :: fcur,fpp
 c
-c     EC power and current densities
+c     effjpl = <J_parallel>/<p_d> /(B_min/<B>)  [A m /W]
+c
+      apar(1) = yg
+      apar(2) = anpl
+      apar(3) = amu
+      apar(4) = anprre
+      apar(5) = dble(ex)
+      apar(6) = dimag(ex)
+      apar(7) = dble(ey)
+      apar(8) = dimag(ey)
+      apar(9) = dble(ez)
+      apar(10) = dimag(ez)
+      apar(11) = dble(ithn)
+c
+      npar=nfpp+necp
+      apar(nfpp+1:npar) = eccdpar(1:necp)
 c
       anpl2=anpl*anpl
-      dnl=1.0_wp_-anpl2
-      ygn=nhn*yg
-      ygn2=ygn*ygn
-
-      resj=0.0_wp_
-      resj1=0.0_wp_
-      resj2=0.0_wp_
-      resp=0.0_wp_
 c
-      rdu2=anpl2+ygn2-1.0_wp_
-      uplp=0.0_wp_
-      uplm=0.0_wp_
-      upltp=0.0_wp_
-      upltm=0.0_wp_
+      effjcd=0.0_wp_
+      anum=0.0_wp_
+      denom=0.0_wp_
+      iokhawa=0
+      ierr=0
+      do nhn=nhmn,nhmx
+        ygn=nhn*yg
+        ygn2=ygn*ygn
 c
-      if (rdu2.ge.0.0_wp_) then
+        rdu2=anpl2+ygn2-1.0_wp_
+c
+        if (rdu2.lt.0.0_wp_) cycle
         rdu=sqrt(rdu2)
+        dnl=1.0_wp_-anpl2
         uplp=(anpl*ygn+rdu)/dnl
         uplm=(anpl*ygn-rdu)/dnl
 c
@@ -5922,92 +5941,72 @@ c
         if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl
         duu=abs(uu1-uu2)
 c
-        apar(1) = yg
-        apar(2) = anpl
-        apar(3) = amu
-        apar(4) = anprre
-        apar(5) = dble(ex)
-        apar(6) = dimag(ex)
-        apar(7) = dble(ey)
-        apar(8) = dimag(ey)
-        apar(9) = dble(ez)
-        apar(10) = dimag(ez)
-        apar(11) = dble(ithn)
+        if(duu.le.dumin) cycle
+c
         apar(12) = dble(nhn)
-        apar(13) = uplp
-        apar(14) = uplm
-        apar(15) = ygn
+        apar(13) = ygn
 c
-        npar=nfpp+necp
-        apar(nfpp+1:npar) = eccdpar(1:necp)
-c
-        if(duu.gt.1.0e-6_wp_) then
-          call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp,
-     .             epp,neval,ier,liw,lw,last,iw,w)
-          if (ier.gt.0) ierr=90
+        call dqagsmv(fpp,uu1,uu2,apar(1:nfpp),nfpp,epsa,epsr,resp,
+     .               epp,neval,ier,liw,lw,last,iw,w)
+        if (ier.gt.0) then
+          ierr=90
+          return
         end if
 c
         rdu2t=cst2*anpl2+ygn2-1.0_wp_
 c
-        if (rdu2t.lt.0.0_wp_.or.cst2.eq.0.0_wp_) then
-c
-c       resonance curve does not cross the trapping region
-c
-          iokhawa=0
-          if(duu.gt.1.0e-4_wp_) then
-            call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
-     .               resj,ej,neval,ier,liw,lw,last,iw,w)
-            if (ier.gt.0) ierr=91
-          end if
-        else
+        if (rdu2t.gt.0.0_wp_.and.cst2.gt.0.0_wp_) then
 c
 c       resonance curve crosses the trapping region
 c
           iokhawa=1
-          rdut=sqrt(rdu2t)
-          cst=sqrt(cst2)
-          upltm=(cst2*anpl*ygn-cst*rdut)/(1.0_wp_-cst2*anpl2)
-          upltp=(cst2*anpl*ygn+cst*rdut)/(1.0_wp_-cst2*anpl2)
+          cstrdut=sqrt(cst2*rdu2t)
+          upltm=(cst2*anpl*ygn-cstrdut)/(1.0_wp_-cst2*anpl2)
+          upltp=(cst2*anpl*ygn+cstrdut)/(1.0_wp_-cst2*anpl2)
+          uleft(0)=uplm
+          uright(0)=upltm
+          uleft(1)=upltp
+          uright(1)=uplp
+        else
 c
-          uu1=uplm
-          uu2=upltm
-          xx1=amu*(anpl*uu1+ygn-1.0_wp_)
-          xx2=amu*(anpl*uu2+ygn-1.0_wp_)
-          if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
-            if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl
-            if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl
-            duu=abs(uu1-uu2)
-            if(duu.gt.1.0e-6_wp_) then
-              call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
-     .        resj1,ej1,neval,ier,liw,lw,last,iw,w)
-              if (ier.gt.0)   then
-              if (abs(resj1).lt.1.0e-10_wp_) then
-                resj1=0.0_wp_
-              else
-                ierr=92
-              end if
-              end if
-            end if
-          end if
+c       resonance curve does not cross the trapping region
 c
-          uu1=upltp
-          uu2=uplp
-          xx1=amu*(anpl*uu1+ygn-1.0_wp_)
-          xx2=amu*(anpl*uu2+ygn-1.0_wp_)
-          if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
-            if(xx2.gt.xxcr) uu2=(xxcr/amu-ygn+1.0_wp_)/anpl
-            if(xx1.gt.xxcr) uu1=(xxcr/amu-ygn+1.0_wp_)/anpl
-            duu=abs(uu1-uu2)
-            if(duu.gt.1.0e-6_wp_) then
-              call dqagsmv(fcur,uu1,uu2,apar,npar,epsa,epsr,
-     .        resj2,ej2,neval,ier,liw,lw,last,iw,w)
-              if (ier.gt.0) then
-                if(ier.ne.2) ierr=93
-            end if
-            end if
-          end if
-          resj=resj1+resj2
+          iokhawa=0
+          uleft(0)=uplm
+          uright(0)=uplp
         end if
+c
+        resj=0.0_wp_
+        do i=0,1
+          xx1=amu*(anpl*uleft(i)+ygn-1.0_wp_)
+          xx2=amu*(anpl*uright(i)+ygn-1.0_wp_)
+          if(xx1.lt.xxcr.or.xx2.lt.xxcr) then
+            if(xx2.gt.xxcr) uright(i)=(xxcr/amu-ygn+1.0_wp_)/anpl
+            if(xx1.gt.xxcr) uleft(i)=(xxcr/amu-ygn+1.0_wp_)/anpl
+            duu=abs(uleft(i)-uright(i))
+            if(duu.gt.dumin) then
+              call dqagsmv(fcur,uleft(i),uright(i),apar,npar,epsa,epsr,
+     .          resji,eji,neval,ier,liw,lw,last,iw,w)
+              if (ier.gt.0)   then
+                if (abs(resji).lt.1.0e-10_wp_) then
+                  resji=0.0_wp_
+                else
+                  ierr=91+iokhawa+i
+                  return
+                end if
+              end if
+            end if
+          end if
+          resj=resj+resji
+          if(iokhawa.eq.0) exit
+        end do
+        anum=anum+resj
+        denom=denom+resp
+      end do
+c
+      if(denom.gt.0.0_wp_) then
+        anucc=canucc*dens*(48.0_wp_-log(1.0e7_wp_*dens*mc2m2*amu**2))
+        effjcd=-ceff*anum/(anucc*denom)
       end if
 
       return
@@ -6035,9 +6034,7 @@ c      extrapar(9) = Re(ez)
 c      extrapar(10) = Im(ez)
 c      extrapar(11) = double(ithn)
 c      extrapar(12) = double(nhn)
-c      extrapar(13) = uplp
-c      extrapar(14) = uplm
-c      extrapar(15) = ygn
+c      extrapar(13) = ygn
 c
       use const_and_precisions, only : wp_,ui=>im
       use math, only : fact
@@ -6048,7 +6045,7 @@ c arguments
       real(wp_), dimension(npar) :: extrapar
 c local variables
       integer :: ithn,nhn,nm,np
-      real(wp_) :: yg,anpl,amu,anprre,uplp,uplm,ygn,upr,upr2,gam,ee,
+      real(wp_) :: yg,anpl,amu,anprre,ygn,upr,upr2,gam,ee,
      . thn2,thn2u,bb,cth,ajbnm,ajbnp,ajbn
       complex(wp_) :: ex,ey,ez,emxy,epxy
 c
@@ -6061,23 +6058,21 @@ c
       ez=cmplx(extrapar(9),extrapar(10),wp_)
       ithn=int(extrapar(11))
       nhn=int(extrapar(12))
-      uplp=extrapar(13)
-      uplm=extrapar(14)
-      ygn=extrapar(15)
+      ygn=extrapar(13)
       
-      upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm)
       gam=anpl*upl+ygn
+      upr2=gam*gam-1.0_wp_-upl*upl
       ee=exp(-amu*(gam-1))
 
-      thn2=1.0_wp_
-      thn2u=upr2*thn2
+!      thn2=1.0_wp_
+      thn2u=upr2 !*thn2
       if(ithn.gt.0) then
         emxy=ex-ui*ey
         epxy=ex+ui*ey
         if(upr2.gt.0.0_wp_) then
           upr=sqrt(upr2)
           bb=anprre*upr/yg
-          if(ithn.eq.1)  then
+          if(ithn.eq.1) then
 c     Larmor radius expansion polarization term at lowest order
             cth=1.0_wp_
             if(nhn.gt.1) cth=(0.5_wp_*bb)**(nhn-1)*nhn/fact(nhn)
@@ -6110,31 +6105,19 @@ c
       function fjch(upl,extrapar,npar)
 c     integration variable upl passed explicitly. Other variables passed
 c     as array of extra parameters of length npar=size(extrapar).
-c     variables with index 1..15 must be passed to fpp
-c     variable with index 16 is zeff
-c     variables with index gt 16 are specific of the cd model
+c     variables with index 1..13 must be passed to fpp
+c     variable with index 14 is zeff
+c     variables with index gt 14 are specific of the cd model
 c     
-c      extrapar(1) = yg
 c      extrapar(2) = anpl
-c      extrapar(3) = amu
 c      extrapar(4) = Re(anprw)
-c      extrapar(5) = Re(ex)
-c      extrapar(6) = Im(ex)
-c      extrapar(7) = Re(ey)
-c      extrapar(8) = Im(ey)
-c      extrapar(9) = Re(ez)
-c      extrapar(10) = Im(ez)
-c      extrapar(11) = double(ithn)
-c      extrapar(12) = double(nhn)
-c      extrapar(13) = uplp
-c      extrapar(14) = uplm
-c      extrapar(15) = ygn
+c      extrapar(13) = ygn
 c
-c      extrapar(16) = zeff
-c      extrapar(17) = rb
-c      extrapar(18) = alams
-c      extrapar(19) = pa
-c      extrapar(20) = fp0s
+c      extrapar(14) = zeff
+c      extrapar(15) = rb
+c      extrapar(16) = alams
+c      extrapar(17) = pa
+c      extrapar(18) = fp0s
 c
       use const_and_precisions, only : wp_
       use conical, only : fconic
@@ -6145,25 +6128,23 @@ c arguments
       real(wp_), dimension(npar) :: extrapar
 c local variables
       integer :: nhn
-      real(wp_) :: anpl,anprre,uplp,uplm,ygn,zeff,rb,alams,pa,fp0s,
+      real(wp_) :: anpl,anprre,ygn,zeff,rb,alams,pa,fp0s,
      . upr2,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,alam,fp0,
      . dfp0,fh,dfhl,eta,fpp
       complex(wp_) :: ex,ey,ez
 c
       anpl=extrapar(2)
       anprre=extrapar(4)
-      uplp=extrapar(13)
-      uplm=extrapar(14)
-      ygn=extrapar(15)
-      zeff=extrapar(16)
-      rb=extrapar(17)
-      alams=extrapar(18)
-      pa=extrapar(19)
-      fp0s=extrapar(20)
+      ygn=extrapar(13)
+      zeff=extrapar(14)
+      rb=extrapar(15)
+      alams=extrapar(16)
+      pa=extrapar(17)
+      fp0s=extrapar(18)
 
-      upr2=(1.0_wp_-anpl**2)*(uplp-upl)*(upl-uplm)
       gam=anpl*upl+ygn
       u2=gam*gam-1.0_wp_
+      upr2=u2-upl*upl
       u=sqrt(u2)
       z5=Zeff+5.0_wp_
       xi=1.0_wp_/z5**2
@@ -6187,7 +6168,7 @@ c
       eta=gam*fh*(gg/u+dgg)+upl*(anpl*u2-upl*gam)*gg*dfhl/(u2*u*rb*alam)
 
       if(upl.lt.0.0_wp_) eta=-eta
-      fjch=eta*fpp(upl,extrapar,npar)
+      fjch=eta*fpp(upl,extrapar(1:13),13)
 c
       return
       end
@@ -6201,23 +6182,10 @@ c     variables with index 1..15 must be passed to fpp
 c     variable with index 16 is zeff
 c     variables with index gt 16 are specific of the cd model
 c     
-c      extrapar(1) = yg
 c      extrapar(2) = anpl
-c      extrapar(3) = amu
-c      extrapar(4) = Re(anprw)
-c      extrapar(5) = Re(ex)
-c      extrapar(6) = Im(ex)
-c      extrapar(7) = Re(ey)
-c      extrapar(8) = Im(ey)
-c      extrapar(9) = Re(ez)
-c      extrapar(10) = Im(ez)
-c      extrapar(11) = double(ithn)
-c      extrapar(12) = double(nhn)
-c      extrapar(13) = uplp
-c      extrapar(14) = uplm
-c      extrapar(15) = ygn
+c      extrapar(13) = ygn
 c
-c      extrapar(16) = zeff
+c      extrapar(14) = zeff
 c
       use const_and_precisions, only : wp_
       implicit none
@@ -6232,8 +6200,8 @@ c local variables
       complex(wp_) :: ex,ey,ez
 c
       anpl=extrapar(2)
-      ygn=extrapar(15)
-      zeff=extrapar(16)
+      ygn=extrapar(13)
+      zeff=extrapar(14)
 
       gam=anpl*upl+ygn
       u2=gam*gam-1.0_wp_
@@ -6248,7 +6216,7 @@ c
       gg=u*gu/z5
       dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
       eta=anpl*gg+gam*upl*dgg/u
-      fjch0=eta*fpp(upl,extrapar,npar)
+      fjch0=eta*fpp(upl,extrapar(1:13),13)
 c
       return
       end
@@ -6258,31 +6226,19 @@ c
       function fjncl(upl,extrapar,npar)
 c     integration variable upl passed explicitly. Other variables passed
 c     as array of extra parameters of length npar=size(extrapar).
-c     variables with index 1..15 must be passed to fpp
-c     variable with index 16 is zeff
-c     variables with index gt 16 are specific of the cd model
+c     variables with index 1..13 must be passed to fpp
+c     variable with index 14 is zeff
+c     variables with index gt 14 are specific of the cd model
 c     
-c      extrapar(1) = yg
 c      extrapar(2) = anpl
 c      extrapar(3) = amu
-c      extrapar(4) = Re(anprw)
-c      extrapar(5) = Re(ex)
-c      extrapar(6) = Im(ex)
-c      extrapar(7) = Re(ey)
-c      extrapar(8) = Im(ey)
-c      extrapar(9) = Re(ez)
-c      extrapar(10) = Im(ez)
-c      extrapar(11) = double(ithn)
-c      extrapar(12) = double(nhn)
-c      extrapar(13) = uplp
-c      extrapar(14) = uplm
-c      extrapar(15) = ygn
+c      extrapar(13) = ygn
 c
-c      extrapar(16) = zeff
-c      extrapar(17) = fc
-c      extrapar(18) = rbx
-c      extrapar(19:18+(npar-18)/2) = tlm
-c      extrapar(19+(npar-18)/2:npar) = chlm
+c      extrapar(14) = zeff
+c      extrapar(15) = fc
+c      extrapar(16) = rbx
+c      extrapar(17:16+(npar-16)/2) = tlm
+c      extrapar(17+(npar-16)/2:npar) = chlm
 c
       use const_and_precisions, only : wp_
       use  green_func_p, only: GenSpitzFunc
@@ -6300,15 +6256,15 @@ c local variables
      . bth,uth,fk,dfk,alam,fu,dfu,eta,fpp
 c local variables
       integer :: ier
-      real(wp_), dimension((npar-18)/2) :: wrk
+      real(wp_), dimension((npar-16)/2) :: wrk
       real(wp_), dimension(1) :: xs,ys
 c
       anpl=extrapar(2)
       amu=extrapar(3)
-      ygn=extrapar(15)
-      zeff=extrapar(16)
-      fc=extrapar(17)
-      rbx=extrapar(18)
+      ygn=extrapar(13)
+      zeff=extrapar(14)
+      fc=extrapar(15)
+      rbx=extrapar(16)
 
       gam=anpl*upl+ygn
       u2=gam*gam-1.0_wp_
@@ -6321,52 +6277,26 @@ c
       dfk=dfk*(2.0_wp_/amu)*bth
       
       alam=upr2/u2/rbx
-      nlm=(npar-18)/2
       xs(1)=alam
+      nlm=(npar-16)/2
 c
-      call splev(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3,
+c      extrapar(17:16+(npar-16)/2) = tlm
+c      extrapar(17+(npar-16)/2:npar) = chlm
+c
+      call splev(extrapar(17:16+nlm),nlm,extrapar(17+nlm:npar),3,
      .           xs(1),ys(1),1,ier)
       fu=ys(1)
-      call splder(extrapar(19:18+nlm),nlm,extrapar(19+nlm:npar),3,1,
+      call splder(extrapar(17:16+nlm),nlm,extrapar(17+nlm:npar),3,1,
      .            xs(1),ys(1),1,wrk,ier)
       dfu=ys(1)
       
       eta=gam*fu*dfk/u-2.0_wp_*(anpl-gam*upl/u2)*fk*dfu*upl/u2/rbx     
       if(upl.lt.0) eta=-eta
-      fjncl=eta*fpp(upl,extrapar,npar)
+      fjncl=eta*fpp(upl,extrapar(1:13),13)
       return
       end
 c
 c
-c
-C      subroutine vlambda(alam,tlm,chlm,nlmt,fv,dfv)
-C      use const_and_precisions, only : wp_
-C      use dierckx, only : splev,splder
-C      implicit none
-Cc local constants
-C      integer, parameter :: ksp=3
-Cc arguments
-C      integer :: nlmt
-C      real(wp_) :: alam,fv,dfv
-C      real(wp_), dimension(nlmt) :: tlm,chlm
-Cc local variables
-C      integer :: nlm,ier
-C      real(wp_), dimension(nlmt) :: wrk
-C      real(wp_), dimension(1) :: xxs,ffs
-Cc
-C      nlm=nlmt
-C      xxs(1)=alam
-Cc
-C      call splev(tlm,nlm,chlm,ksp,xxs(1),ffs(1),1,ier)
-C      fv=ffs(1)
-Cc
-C      call splder(tlm,nlm,chlm,ksp,1,xxs(1),ffs(1),1,wrk,ier)
-C      dfv=ffs(1)
-Cc
-C      return
-C      end  
-c
-c
 c
       subroutine projxyzt(iproj,nfile)
       use const_and_precisions, only : wp_