refl/gray/after_onestep: fixed test for search of ray that completed 1st pass in plasma

Makefile

@@ -2,40 +2,63 @@
 EXE=gray
 
 # Objects list
-OBJ=gray.o grayl.o reflections.o green_func_p.o \
-    const_and_precisions.o itm_constants.o itm_types.o
-
+MAINOBJ=gray.o
+OTHOBJ=conical.o const_and_precisions.o coreprofiles.o dierckx.o dispersion.o \
+     eccd.o eierf.o graydata_anequil.o graydata_flags.o graydata_par.o \
+    equilibrium.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
+    reflections.o simplespline.o utils.o beamdata.o 
+    
 # Alternative search paths
 vpath %.f90 src
 vpath %.f src
 
 # Fortran compiler name and flags
 FC=gfortran
-FFLAGS=-Wall -fcheck=all
+FFLAGS=-Wall -g -fcheck=all
 
 DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
 
 all: $(EXE)
 
 # Build executable from object files
-$(EXE): $(OBJ)
+$(EXE): $(MAINOBJ) $(OTHOBJ)
 	$(FC) $(FFLAGS) -o $@ $^
 
 # Dependencies on modules
-gray.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
+gray.o: const_and_precisions.o coreprofiles.o dierckx.o dispersion.o eccd.o \
+        graydata_anequil.o graydata_flags.o graydata_par.o \
+        equilibrium.o magsurf_data.o math.o numint.o quadpack.o \
+        reflections.o simplespline.o utils.o beamdata.o 
+conical.o: const_and_precisions.o
+coreprofiles.o: const_and_precisions.o dierckx.o graydata_anequil.o \
+                graydata_flags.o simplespline.o utils.o
+dierckx.o: const_and_precisions.o
+dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
+eccd.o: const_and_precisions.o conical.o magsurf_data.o dierckx.o numint.o
+eierf.o: const_and_precisions.o
+graydata_anequil.o: const_and_precisions.o
+graydata_flags.o: const_and_precisions.o
+graydata_par.o: const_and_precisions.o
+equilibrium.o: const_and_precisions.o dierckx.o minpack.o simplespline.o utils.o
+magsurf_data.o: const_and_precisions.o
+math.o: const_and_precisions.o
+minpack.o: const_and_precisions.o
+numint.o: const_and_precisions.o
+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
 	$(FC) $(FFLAGS) -c $<
 
-gray.o:gray.f green_func_p.o
-	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
+%.o: %.f
+	$(FC) $(FFLAGS) -c $<
 
-grayl.o:grayl.f
-	$(FC) $(FFLAGS) -c $^ 
+gray.o:gray.f
+	$(FC) -cpp $(DIRECTIVES) $(FFLAGS) -c $< 
 
 .PHONY: clean install
 # Remove output files

input/gray.data

@@ -0,0 +1,28 @@
+40. 0.           : alfac, betac poloidal and toroidal launching angles 
+82.7              : f (GHz)
+1                 : P0 (MW)   injected power
+1 1 1             : nray ktx,rhomx
+123. 0. 50.       : x00  y00 z00 (cm)  mirror position (y00=0 !)
+2.1 2.1 162 162 0 : w0xt w0yt pw0xt pw0yt awr  (cm)  [pw0>0 towards plasma]
+0                : ibeam=0 uses data above, >0 read from file beam data (1,2,3)
+mpusm            : filename.txt for beam data max(character*16)
+2     0          : iequil= 0 vacuum, 1 analytical equilibrium, 2 EQDSK;indXpoint
+0                : iprof= 0/1 analytical/numerical dens & temp profiles
+1 2 -10           : iwarm= 0-3 ECRH&CD computation,  ilarm  order of larmor expans
+0                : ieccd 0/1 NO/YES  ECCD calculation
+1  251          : ipec=0/1 profiles in dpsi/drhop; nnd=no. intervals +1 (=< 501)
+.01 50000 100 5 1  : dst (cm), nsteps < = 8000, istprj,istpl,ist (0,1,2 integration in s,ct,Sr)
+0 -2 0.624       : igrad (0  rtr - 1 gauss beam), ipass=1/2 passes in plasma, R_wall(m)
+2 0 0 0         : IOX = 1  :  O-MODE  IOX = 2  :  X-MODE, psi_pol_ant chi_pol_ant
+TCV             : filename EQDSK : filename.eqdsk max(character 16)
+0 0.001 0 1     : ipsinorm =0/1 [psi normalized (1) or not (0) in EQDSK],  sspl spline coeff psi
+1 1 1 1         : factb,  iscal=0/1/2 nustar=const/ngreenw=const/no_scal, factT factn scal.
+filename        : filename profiles : filename.prf  max(character*16)
+1.015           : psi plasma boundary
+-1 -1 2        : signum B_phi I_phi +1=counterclockwise , -1=clockwise from above, icocos (0 old case, 11 iter)
+9.99999999999999999992 8 8   : dens190 aln1,aln2   (10^19 m-3)
+4.99999999999999999424 0 8 8     : te0,dte0,alt1,alt2  (keV)
+2               : Zeff
+620 0 180       : rr0 zr0 a (cm)
+5.3             : b0 (T)
+1 3  2          : q0,qa alq

src/beamdata.f90

@@ -0,0 +1,80 @@
+module beamdata
+  use const_and_precisions, only : wp_
+  implicit none
+
+  integer, parameter :: jmx=31,kmx=36,nmx=100000
+  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

src/conical.f90

@@ -0,0 +1,853 @@
+module conical
+
+  use const_and_precisions, only : wp_
+  implicit none
+
+contains
+
+  function fconic(x,tau,m)
+!
+! this function subprogram computes the conical functions of the
+! first kind  P sub(-1/2 + i*tau) (x)  for m = 0 and m = 1.
+! Ref. in Kolbig, Comp. Phys. Commun. 23 (1981) 51
+!
+    implicit none
+    real(wp_), intent(in) :: x, tau
+    integer, intent(in) :: m
+    real(wp_) :: fconic
+    real(wp_) :: t(7),h(9),v(11)
+    real(wp_) :: aa,a0,a1,a2,a3,a4,a5,a6,b0,b1,fm,fn,fn1,r1,r2,s,s0,s1
+    real(wp_) :: x1,y,y2,y3,z
+    integer :: jp,j,n
+    real(wp_), parameter :: rpi=1.7724538509055_wp_,pi2=0.63661977236758_wp_
+    real(wp_), parameter :: eps=1.0e-14_wp_
+    integer, parameter :: nout=2,nmax=200
+!
+    complex(wp_) a,b,c,ti,r,rr,q,u,u0,u1,u2,uu
+    complex(wp_) v0,v1,v2,vv,w(19)
+!
+    logical lm0,lm1,lta
+
+    fconic=0.0_wp_
+    lm0=m == 0
+    lm1=m == 1
+    if(.not.(lm0 .or. lm1)) then
+      write(nout,"(1x,'fconic ... illegal value for m = ',i4)") m
+      return
+    end if
+    fm=m
+    fconic=1.0_wp_-fm
+    if(x == 1.0_wp_) return
+!
+    fconic=0.0_wp_
+    if(tau == 0.0_wp_ .and. abs(x-1.0_wp_) > 0.01_wp_) then
+      if(x > 1.0_wp_) then
+        y=sqrt((x-1.0_wp_)/(x+1.0_wp_))
+        z=ellick(y)
+        s=sqrt(0.5_wp_*(x+1.0_wp_))
+        if(lm0) fconic=pi2*z/s
+        if(lm1) fconic=pi2*s*(ellice(y)-z)/sqrt(x**2-1.0_wp_)
+        return
+      else
+        y=sqrt(0.5_wp_*(1.0_wp_-x))
+        z=ellick(y)
+        if(lm0) fconic=pi2*z
+        if(lm1) fconic=pi2*(ellice(y)-0.5_wp_*(1.0_wp_+x)*z)/ &
+                       sqrt(1.0_wp_-x**2)
+        return
+      end if
+    else
+      ti=cmplx(0._wp_,tau,wp_)
+!
+      if((-1._wp_ < x .and. x <= 0.0_wp_).or. &
+         (0.0_wp_ < x .and. x <= 0.1_wp_ .and.tau<= 17.0_wp_).or. &
+         (0.1_wp_ < x .and. x <= 0.2_wp_ .and.tau<=  5.0_wp_)) then
+        lta=tau <= 10.0_wp_
+        x1=x**2
+        a=0.5_wp_*(0.5_wp_-fm-ti)
+        b=0.5_wp_*(0.5_wp_-fm+ti)
+        c=0.5_wp_
+        jp=30
+      else if((0.1_wp_ < x .and. x <= 0.2_wp_ .and.tau<= 17.0_wp_) &
+          .or.(0.2_wp_ < x .and. x <= 1.5_wp_ .and.tau<= 20.0_wp_)) &
+          then
+        lta=x > 1.0_wp_ .or. x <= 1.0_wp_ .and. tau <= 5.0_wp_
+        x1=(1.0_wp_-x)/2._wp_
+        a=0.5_wp_+fm-ti
+        b=0.5_wp_+fm+ti
+        c=fm+1.0_wp_
+        jp=32
+      else if(1.5_wp_ < x .and. tau <= max(20.0_wp_,x)) then
+        lta=.true.
+        x1=1.0_wp_/x**2
+        u=exp((-0.5_wp_+ti)*log(2.0_wp_*x)+clogam(1.0_wp_+ti) &
+          -clogam(1.5_wp_-fm+ti))
+        a=0.5_wp_*(0.5_wp_-fm-ti)
+        b=0.5_wp_*(1.5_wp_-fm-ti)
+        c=1.0_wp_-ti
+        jp=33
+      else
+        if(x > 1.0_wp_) then
+          s=sqrt(x**2-1.0_wp_)
+          t(1)=log(x+s)
+          h(1)=tau*t(1)
+          b0=besj0l(h(1))
+          b1=besj1l(h(1))
+          z=1.0_wp_
+        else
+          s=sqrt(1.0_wp_-x**2)
+          t(1)=acos(x)
+          h(1)=tau*t(1)
+          b0=besi0(h(1))
+          b1=besi1(h(1))
+          z=-1.0_wp_
+        end if
+        h(1)=t(1)*x/s
+        v(1)=tau
+        do j = 2,7
+          t(j)=t(j-1)*t(1)
+          h(j)=h(j-1)*h(1)
+        end do
+        do j = 2,11
+          v(j)=v(j-1)*v(1)
+        end do
+!
+        if(lm1) then
+          aa=-1.0_wp_
+          a0=3.0_wp_*(1.0_wp_-h(1))/(8.0_wp_*t(1))
+          a1=(-15.0_wp_*h(2)+6.0_wp_*h(1)+9.0_wp_+z*8.0_wp_*t(2))/ &
+              (128.0_wp_*t(2))
+          a2=3.0_wp_*(-35.0_wp_*h(3)-15.0_wp_*h(2)+15.0_wp_*h(1)+35.0_wp_ &
+             +z*t(2)*(32.0_wp_*h(1)+8.0_wp_))/(1024.0_wp_*t(3))
+          a3=(-4725.0_wp_*h(4)-6300.0_wp_*h(3)-3150.0_wp_*h(2)+3780.0_wp_*h(1) &
+              +10395.0_wp_-1216.0_wp_*t(4)+z*t(2)*(6000.0_wp_*h(2) &
+              +5760.0_wp_*h(1)+1680.0_wp_)) /(32768.0_wp_*t(4))
+          a4=7.0_wp_*(-10395.0_wp_*h(5)-23625.0_wp_*h(4)-28350.0_wp_*h(3) &
+             -14850.0_wp_*h(2)+19305.0_wp_*h(1)+57915.0_wp_ &
+             -t(4)*(6336.0_wp_*h(1)+6080.0_wp_)+z*t(2)*(16800.0_wp_*h(3) &
+             +30000.0_wp_*h(2)+25920.0_wp_*h(1)+7920.0_wp_))/ &
+             (262144.0_wp_*t(5))
+          a5=(-2837835.0_wp_*h(6)-9168390.0_wp_*h(5)-16372125.0_wp_*h(4) &
+              -18918900*h(3) -10135125.0_wp_*h(2)+13783770.0_wp_*h(1) &
+              +43648605.0_wp_-t(4)*(3044160.0_wp_*h(2)+5588352.0_wp_*h(1) &
+              +4213440.0_wp_)+z*t(2)*(5556600.0_wp_*h(4)+14817600.0_wp_*h(3) &
+              +20790000.0_wp_*h(2)+17297280.0_wp_*h(1)+5405400.0_wp_ &
+              +323072.0_wp_*t(4)))/ (4194304.0_wp_*t(6))
+          a6=0.0_wp_
+        else
+          aa=0.0_wp_
+          a0=1.0_wp_
+          a1=(h(1)-1.0_wp_)/(8.0_wp_*t(1))
+          a2=(9.0_wp_*h(2)+6.0_wp_*h(1)-15.0_wp_-z*8.0_wp_*t(2))/ &
+              (128.0_wp_*t(2))
+          a3=5.0_wp_*(15.0_wp_*h(3)+27.0_wp_*h(2)+21.0_wp_*h(1)-63.0_wp_ &
+             -z*t(2)*(16.0_wp_*h(1)+24.0_wp_))/(1024.0_wp_*t(3))
+          a4=7.0_wp_*(525.0_wp_*h(4)+1500.0_wp_*h(3)+2430.0_wp_*h(2) &
+             +1980.0_wp_*h(1)-6435.0_wp_+192.0_wp_*t(4)-z*t(2)* &
+              (720.0_wp_*h(2)+1600.0_wp_*h(1)+2160.0_wp_))/(32768.0_wp_*t(4))
+          a5=21.0_wp_*(2835.0_wp_*h(5)+11025.0_wp_*h(4)+24750.0_wp_*h(3) &
+             +38610.0_wp_*h(2)+32175.0_wp_*h(1)-109395.0_wp_+t(4) &
+             *(1984.0_wp_*h(1)+4032.0_wp_)-z*t(2) &
+             *(4800.0_wp_*h(3)+15120.0_wp_*h(2)+26400.0_wp_*h(1)+34320.0_wp_)) &
+             /(262144.0_wp_*t(5))
+          a6=11.0_wp_*(218295.0_wp_*h(6)+1071630.0_wp_*h(5)+3009825.0_wp_*h(4) &
+             +6142500.0_wp_*h(3)+9398025.0_wp_*h(2)+7936110.0_wp_*h(1) &
+             -27776385.0_wp_+t(4)*(254016.0_wp_*h(2) &
+             +749952.0_wp_*h(1)+1100736.0_wp_)-z*t(2)*(441000.0_wp_*h(4) &
+             +1814400.0_wp_*h(3)+4127760.0_wp_*h(2)+6552000.0_wp_*h(1) &
+             +8353800.0_wp_+31232.0_wp_*t(4)))/(4194304.0_wp_*t(6))
+        end if
+        s0=a0+(-4.0_wp_*a3/t(1)+a4)/v(4)+(-192.0_wp_*a5/t(3) &
+           +144.0_wp_*a6/t(2))/v(8)+z*(-a2/v(2)+(-24.0_wp_*a4/t(2) &
+           +12.0_wp_*a5/t(1)-a6)/v(6)+(-1920.0_wp_*a6/t(4))/v(10))
+        s1=a1/v(1)+(8.0_wp_*(a3/t(2)-a4/t(1))+a5)/v(5)+(384.0_wp_*a5/t(4) &
+           -768.0_wp_*a6/t(3))/v(9)+z*(aa*v(1)+(2.0_wp_*a2/t(1)-a3)/v(3) &
+           +(48.0_wp_*a4/t(3)-72.0_wp_*a5/t(2) &
+           +18.0_wp_*a6/t(1))/v(7)+(3840.0_wp_*a6/t(5))/v(11))
+        fconic=sqrt(t(1)/s)*(b0*s0+b1*s1)
+        return
+      end if
+!
+      do
+        if(lta) then
+          y=-x1
+          y2=y**2
+          y3=y**3
+          w(1)=a+1.0_wp_
+          w(2)=a+2.0_wp_
+          w(3)=b+1.0_wp_
+          w(4)=b+2.0_wp_
+          w(5)=c+1.0_wp_
+          w(6)=c*w(5)
+          w(7)=a+b
+          w(8)=a*b
+          w(9)=(w(8)/c)*y
+          w(10)=w(1)*w(3)
+          w(11)=w(2)*w(4)
+          w(12)=1.0_wp_+(w(11)/(2.0_wp_*w(5)))*y
+          w(13)=w(7)-6.0_wp_
+          w(14)=w(7)+6.0_wp_
+          w(15)=2.0_wp_-w(8)
+          w(16)=w(15)-2.0_wp_*w(7)
+!
+          v0=1.0_wp_
+          v1=1.0_wp_+(w(10)/(2.0_wp_*c))*y
+          v2=w(12)+(w(10)*w(11)/(12.0_wp_*w(6)))*y2
+          u0=1.0_wp_
+          u1=v1-w(9)
+          u2=v2-w(9)*w(12)+(w(8)*w(10)/(2.0_wp_*w(6)))*y2
+!
+          r=1.0_wp_
+          n=2
+          do
+            n=n+1
+            if(n > nmax) then
+              write(nout,200) x,tau,m
+              return
+            end if
+            rr=r
+            fn=n
+            h(1)=fn-1.0_wp_
+            h(2)=fn-2.0_wp_
+            h(3)=fn-3.0_wp_
+            h(4)=2.0_wp_*fn
+            h(5)=h(4)-3.0_wp_
+            h(6)=2.0_wp_*h(5)
+            h(7)=4.0_wp_*(h(4)-1.0_wp_)*h(5)
+            h(8)=8.0_wp_*h(5)**2*(h(4)-5.0_wp_)
+            h(9)=3.0_wp_*fn**2
+            w(1)=a+h(1)
+            w(2)=a+h(2)
+            w(3)=b+h(1)
+            w(4)=b+h(2)
+            w(5)=c+h(1)
+            w(6)=c+h(2)
+            w(7)=c+h(3)
+            w(8)=h(2)-a
+            w(9)=h(2)-b
+            w(10)=h(1)-c
+            w(11)=w(1)*w(3)
+            w(12)=w(5)*w(6)
+!
+            w(17)=1.0_wp_+((h(9)+w(13)*fn+w(16))/(h(6)*w(5)))*y
+            w(18)=-((w(11)*w(10)/h(6)+(h(9)-w(14)*fn+w(15))* &
+                     w(11)*y/h(7))/w(12))*y
+            w(19)=(w(2)*w(11)*w(4)*w(8)*w(9)/(h(8)*w(7)*w(12)))*y3
+            vv=w(17)*v2+w(18)*v1+w(19)*v0
+            uu=w(17)*u2+w(18)*u1+w(19)*u0
+            r=uu/vv
+            if(abs(r-rr) < eps) exit
+            v0=v1
+            v1=v2
+            v2=vv
+            u0=u1
+            u1=u2
+            u2=uu
+          end do
+        else
+          r=1.0_wp_
+          q=1.0_wp_
+          do n = 1,nmax
+            fn=n
+            fn1=fn-1.0_wp_
+            rr=r
+            q=q*x1*(a+fn1)*(b+fn1)/((c+fn1)*fn)
+            r=r+q
+            if(abs(r-rr) < eps) exit
+          end do
+          if (n > nmax) then
+            write(nout,200) x,tau,m
+            return
+          end if
+        end if
+        if (jp/=30) exit
+        r1=real(r)/abs(exp(clogam(a+0.5_wp_)))**2
+        a=0.5_wp_*(1.5_wp_-fm-ti)
+        b=0.5_wp_*(1.5_wp_-fm+ti)
+        c=1.5_wp_
+        jp=31
+      end do
+      if (jp==31) then
+        r2=real(r)/abs(exp(clogam(a-0.5_wp_)))**2
+        fconic=rpi*(r1-2.0_wp_*x*r2)
+        if(lm1) fconic=(2.0_wp_/sqrt(1.0_wp_-x1))*fconic
+        return
+      else if (jp==32) then
+        fconic=real(r)
+        if(.not.lm0) then
+          fconic=0.5_wp_*(tau**2+0.25_wp_)*sqrt(abs(x**2-1.0_wp_))*fconic
+          if(x > 1.0_wp_) fconic=-fconic
+        end if
+        return
+      else if (jp==33) then
+        fconic=2.0_wp_*real(u*r*(0.5_wp_-fm+ti)/ti)/rpi
+        if(lm1) fconic=fconic/sqrt(1.0_wp_-x1)
+        return
+      end if
+    end if
+!
+    200 format(1x,'fconic ... convergence difficulties for c function, x = ', &
+     e12.4,5x,'tau = ',e12.4,5x,'m = ',i5)
+!
+  end function fconic
+
+  function clogam(z)
+!
+    implicit none
+    complex(wp_) :: clogam
+    complex(wp_), intent(in) :: z
+    complex(wp_) :: v,h,r
+    integer :: i,n
+    real(wp_) :: x,t,a,c,d,e,f
+    integer, parameter :: nout=2
+    real(wp_), parameter :: pi=3.1415926535898_wp_
+    real(wp_), dimension(10), parameter :: b= &
+           (/+8.3333333333333e-2_wp_, -2.7777777777778e-3_wp_, &
+            +7.9365079365079e-4_wp_, -5.9523809523810e-4_wp_, &
+            +8.4175084175084e-4_wp_, -1.9175269175269e-3_wp_, &
+            +6.4102564102564e-3_wp_, -2.9550653594771e-2_wp_, &
+            +1.7964437236883e-1_wp_, -1.3924322169059e+0_wp_/)
+!
+    x=real(z)
+    t=aimag(z)
+    if(-abs(x) == aint(x) .and. t == 0.0_wp_) then
+      write(nout,'(1x,f20.2)') x
+      clogam=(0.0_wp_,0.0_wp_)
+      return
+    end if
+    f=abs(t)
+    v=cmplx(x,f,wp_)
+    if(x < 0.0_wp_) v=1.0_wp_-v
+    h=(0.0_wp_,0.0_wp_)
+    c=real(v)
+    if(c < 7.0_wp_) then
+      n=6-int(c)
+      h=v
+      d=aimag(v)
+      a=atan2(d,c)
+      do i = 1,n
+        c=c+1.0_wp_
+        v=cmplx(c,d,wp_)
+        h=h*v
+        a=a+atan2(d,c)
+      end do
+      h=cmplx(0.5_wp_*log(real(h)**2+aimag(h)**2),a,wp_)
+      v=v+1.0_wp_
+    end if
+    r=1.0_wp_/v**2
+    clogam=0.91893853320467_wp_+(v-0.5_wp_)*log(v)-v+(b(1)+r*(b(2)+r*(b(3) &
+     +r*(b(4)+r*(b(5)+r*(b(6)+r*(b(7)+r*(b(8)+r*(b(9)+r*b(10)))))))))) &
+     /v-h
+    if(x < 0.0_wp_) then
+!
+      a=aint(x)-1.0_wp_
+      c=pi*(x-a)
+      d=pi*f
+      e=exp(-2.0_wp_*d)
+      f=sin(c)
+      e=d+0.5_wp_*log(e*f**2+0.25_wp_*(1.0_wp_-e)**2)
+      f=atan2(cos(c)*tanh(d),f)-a*pi
+      clogam=1.1447298858494_wp_-cmplx(e,f,wp_)-clogam
+!
+    end if
+    if(t < 0.0_wp_) clogam=conjg(clogam)
+!
+  end function clogam
+
+  function ellick(xk)
+    implicit none
+    real(wp_), intent(in) :: xk
+    real(wp_) :: ellick, ellice
+    integer :: i
+    real(wp_) :: eta,pa,pb,pc,pd
+    real(wp_), dimension(10), parameter :: &
+           a=(/9.6573590280856e-2_wp_, 3.0885146271305e-2_wp_, &
+               1.4938013532687e-2_wp_, 8.7898018745551e-3_wp_, &
+               6.1796274460533e-3_wp_, 6.8479092826245e-3_wp_, &
+               9.8489293221769e-3_wp_, 8.0030039806500e-3_wp_, &
+               2.2966348983970e-3_wp_, 1.3930878570066e-4_wp_/), &
+           b=(/1.2499999999991e-1_wp_, 7.0312499739038e-2_wp_, &
+               4.8828041906862e-2_wp_, 3.7377739758624e-2_wp_, &
+               3.0124849012899e-2_wp_, 2.3931913323111e-2_wp_, &
+               1.5530941631977e-2_wp_, 5.9739042991554e-3_wp_, &
+               9.2155463496325e-4_wp_, 2.9700280966556e-5_wp_/), &
+           c=(/4.4314718056089e-1_wp_, 5.6805194567559e-2_wp_, &
+               2.1831811676130e-2_wp_, 1.1569595745295e-2_wp_, &
+               7.5950934225594e-3_wp_, 7.8204040609596e-3_wp_, &
+               1.0770635039866e-2_wp_, 8.6384421736041e-3_wp_, &
+               2.4685033304607e-3_wp_, 1.4946621757181e-4_wp_/), &
+           d=(/2.4999999999990e-1_wp_, 9.3749999721203e-2_wp_, &
+               5.8593661255531e-2_wp_, 4.2717890547383e-2_wp_, &
+               3.3478943665762e-2_wp_, 2.6145014700314e-2_wp_, &
+               1.6804023346363e-2_wp_, 6.4321465864383e-3_wp_, &
+               9.8983328462254e-4_wp_, 3.1859195655502e-5_wp_/)
+!
+    if(abs(xk) >= 1.0_wp_) then
+      ellick=0.0_wp_
+      return
+    end if
+    eta=1.0_wp_-xk**2
+    pa=a(10)
+    do i = 1,9
+      pa=pa*eta+a(10-i)
+    end do
+    pa=pa*eta
+    pb=b(10)
+    do i = 1,9
+      pb=pb*eta+b(10-i)
+    end do
+    pb=pb*eta
+    ellick=1.3862943611199_wp_+pa-log(eta)*(0.5_wp_+pb)
+    return
+!
+    entry ellice(xk)
+!
+    if (abs(xk) >= 1.0_wp_) then
+      if (abs(xk) > 1.0_wp_) then
+        ellick=0.0_wp_
+      else
+        ellick=1.0_wp_
+      end if
+      return
+    end if
+    eta=1.0_wp_-xk**2
+    pc=c(10)
+    do i = 1,9
+      pc=pc*eta+c(10-i)
+    end do
+    pc=pc*eta
+    pd=d(10)
+    do i = 1,9
+      pd=pd*eta+d(10-i)
+    end do
+    pd=pd*eta
+    ellick=1.0_wp_+pc-log(eta)*pd
+  end function ellick
+
+  function besjy(x)
+    implicit none
+    real(wp_), intent(in) :: x
+    real(wp_) :: besjy,besj0l,besj1l
+    real(wp_) :: besy0,besy1
+    logical :: l
+    real(wp_) :: v,f,a,b,p,q
+    integer, parameter :: nout=2
+!
+    entry besj0l(x)
+!
+    l=.true.
+    v=abs(x)
+    if(v >= 8.0_wp_) go to 4
+      8 f=0.0625_wp_*x**2-2.0_wp_
+    a =           - 0.0000000000000008_wp_
+    b = f * a     + 0.0000000000000413_wp_
+    a = f * b - a - 0.0000000000019438_wp_
+    b = f * a - b + 0.0000000000784870_wp_
+    a = f * b - a - 0.0000000026792535_wp_
+    b = f * a - b + 0.0000000760816359_wp_
+    a = f * b - a - 0.0000017619469078_wp_
+    b = f * a - b + 0.0000324603288210_wp_
+    a = f * b - a - 0.0004606261662063_wp_
+    b = f * a - b + 0.0048191800694676_wp_
+    a = f * b - a - 0.0348937694114089_wp_
+    b = f * a - b + 0.1580671023320973_wp_
+    a = f * b - a - 0.3700949938726498_wp_
+    b = f * a - b + 0.2651786132033368_wp_
+    a = f * b - a - 0.0087234423528522_wp_
+    a = f * a - b + 0.3154559429497802_wp_
+    besjy=0.5_wp_*(a-b)
+    if(l) return
+!
+    a =           + 0.0000000000000016_wp_
+    b = f * a     - 0.0000000000000875_wp_
+    a = f * b - a + 0.0000000000040263_wp_
+    b = f * a - b - 0.0000000001583755_wp_
+    a = f * b - a + 0.0000000052487948_wp_
+    b = f * a - b - 0.0000001440723327_wp_
+    a = f * b - a + 0.0000032065325377_wp_
+    b = f * a - b - 0.0000563207914106_wp_
+    a = f * b - a + 0.0007531135932578_wp_
+    b = f * a - b - 0.0072879624795521_wp_
+    a = f * b - a + 0.0471966895957634_wp_
+    b = f * a - b - 0.1773020127811436_wp_
+    a = f * b - a + 0.2615673462550466_wp_
+    b = f * a - b + 0.1790343140771827_wp_
+    a = f * b - a - 0.2744743055297453_wp_
+    a = f * a - b - 0.0662922264065699_wp_
+    besjy=0.636619772367581_wp_*log(x)*besjy+0.5_wp_*(a-b)
+    return
+!
+      4 f=256.0_wp_/x**2-2.0_wp_
+    b =           + 0.0000000000000007_wp_
+    a = f * b     - 0.0000000000000051_wp_
+    b = f * a - b + 0.0000000000000433_wp_
+    a = f * b - a - 0.0000000000004305_wp_
+    b = f * a - b + 0.0000000000051683_wp_
+    a = f * b - a - 0.0000000000786409_wp_
+    b = f * a - b + 0.0000000016306465_wp_
+    a = f * b - a - 0.0000000517059454_wp_
+    b = f * a - b + 0.0000030751847875_wp_
+    a = f * b - a - 0.0005365220468132_wp_
+    a = f * a - b + 1.9989206986950373_wp_
+    p=a-b
+    b =           - 0.0000000000000006_wp_
+    a = f * b     + 0.0000000000000043_wp_
+    b = f * a - b - 0.0000000000000334_wp_
+    a = f * b - a + 0.0000000000003006_wp_
+    b = f * a - b - 0.0000000000032067_wp_
+    a = f * b - a + 0.0000000000422012_wp_
+    b = f * a - b - 0.0000000007271916_wp_
+    a = f * b - a + 0.0000000179724572_wp_
+    b = f * a - b - 0.0000007414498411_wp_
+    a = f * b - a + 0.0000683851994261_wp_
+    a = f * a - b - 0.0311117092106740_wp_
+    q=8.0_wp_*(a-b)/v
+    f=v-0.785398163397448_wp_
+    a=cos(f)
+    b=sin(f)
+    f=0.398942280401432_wp_/sqrt(v)
+    if(l) go to 6
+    besjy=f*(q*a+p*b)
+    return
+      6 besjy=f*(p*a-q*b)
+    return
+!
+    entry besj1l(x)
+!
+    l=.true.
+    v=abs(x)
+    if(v >= 8.0_wp_) go to 5
+      3 f=0.0625_wp_*x**2-2.0_wp_
+    b =           + 0.0000000000000114_wp_
+    a = f * b     - 0.0000000000005777_wp_
+    b = f * a - b + 0.0000000000252812_wp_
+    a = f * b - a - 0.0000000009424213_wp_
+    b = f * a - b + 0.0000000294970701_wp_
+    a = f * b - a - 0.0000007617587805_wp_
+    b = f * a - b + 0.0000158870192399_wp_
+    a = f * b - a - 0.0002604443893486_wp_
+    b = f * a - b + 0.0032402701826839_wp_
+    a = f * b - a - 0.0291755248061542_wp_
+    b = f * a - b + 0.1777091172397283_wp_
+    a = f * b - a - 0.6614439341345433_wp_
+    b = f * a - b + 1.2879940988576776_wp_
+    a = f * b - a - 1.1918011605412169_wp_
+    a = f * a - b + 1.2967175412105298_wp_
+    besjy=0.0625_wp_*(a-b)*x
+    if(l) return
+!
+    b =           - 0.0000000000000244_wp_
+    a = f * b     + 0.0000000000012114_wp_
+    b = f * a - b - 0.0000000000517212_wp_
+    a = f * b - a + 0.0000000018754703_wp_
+    b = f * a - b - 0.0000000568844004_wp_
+    a = f * b - a + 0.0000014166243645_wp_
+    b = f * a - b - 0.0000283046401495_wp_
+    a = f * b - a + 0.0004404786298671_wp_
+    b = f * a - b - 0.0051316411610611_wp_
+    a = f * b - a + 0.0423191803533369_wp_
+    b = f * a - b - 0.2266249915567549_wp_
+    a = f * b - a + 0.6756157807721877_wp_
+    b = f * a - b - 0.7672963628866459_wp_
+    a = f * b - a - 0.1286973843813500_wp_
+    a = f * a - b + 0.0406082117718685_wp_
+    besjy=0.636619772367581_wp_*log(x)*besjy-0.636619772367581_wp_/x &
+         +0.0625_wp_*(a-b)*x
+    return
+!
+      5 f=256.0_wp_/x**2-2.0_wp_
+    b =           - 0.0000000000000007_wp_
+    a = f * b     + 0.0000000000000055_wp_
+    b = f * a - b - 0.0000000000000468_wp_
+    a = f * b - a + 0.0000000000004699_wp_
+    b = f * a - b - 0.0000000000057049_wp_
+    a = f * b - a + 0.0000000000881690_wp_
+    b = f * a - b - 0.0000000018718907_wp_
+    a = f * b - a + 0.0000000617763396_wp_
+    b = f * a - b - 0.0000039872843005_wp_
+    a = f * b - a + 0.0008989898330859_wp_
+    a = f * a - b + 2.0018060817200274_wp_
+    p=a-b
+    b =           + 0.0000000000000007_wp_
+    a = f * b     - 0.0000000000000046_wp_
+    b = f * a - b + 0.0000000000000360_wp_
+    a = f * b - a - 0.0000000000003264_wp_
+    b = f * a - b + 0.0000000000035152_wp_
+    a = f * b - a - 0.0000000000468636_wp_
+    b = f * a - b + 0.0000000008229193_wp_
+    a = f * b - a - 0.0000000209597814_wp_
+    b = f * a - b + 0.0000009138615258_wp_
+    a = f * b - a - 0.0000962772354916_wp_
+    a = f * a - b + 0.0935555741390707_wp_
+    q=8.0_wp_*(a-b)/v
+    f=v-2.356194490192345_wp_
+    a=cos(f)
+    b=sin(f)
+    f=0.398942280401432_wp_/sqrt(v)
+    if(l) go to 7
+    besjy=f*(q*a+p*b)
+    return
+      7 besjy=f*(p*a-q*b)
+    if(x < 0.0_wp_) besjy=-besjy
+    return
+!
+    entry besy0(x)
+!
+    if(x <= 0.0_wp_) go to 9
+    l=.false.
+    v=x
+    if(v >= 8.0_wp_) go to 4
+    go to 8
+    entry besy1(x)
+!
+    if(x <= 0.0_wp_) go to 9
+    l=.false.
+    v=x
+    if(v >= 8.0_wp_) go to 5
+    go to 3
+!
+      9 besjy=0.0_wp_
+    write(nout,"(1x,'besjy ... non-positive argument x = ',e15.4)") x
+  end function besjy
+
+  function besik(x)
+    implicit none
+    real(wp_), intent(in) :: x
+    real(wp_) :: besik,ebesi0,besi0,ebesi1,besi1,ebesk0,besk0,ebesk1,besk1
+    logical :: l,e
+    real(wp_) :: v,f,a,b,z
+    integer, parameter :: nout=2
+!
+    entry ebesi0(x)
+!
+    e=.true.
+    go to 1
+!
+    entry besi0(x)
+!
+    e=.false.
+      1 l=.true.
+    v=abs(x)
+    if(v >= 8.0_wp_) go to 4
+      8 f=0.0625_wp_*x**2-2.0_wp_
+    a = 0.000000000000002_wp_
+    b = f * a     +   0.000000000000120_wp_
+    a = f * b - a +   0.000000000006097_wp_
+    b = f * a - b +   0.000000000268828_wp_
+    a = f * b - a +   0.000000010169727_wp_
+    b = f * a - b +   0.000000326091051_wp_
+    a = f * b - a +   0.000008738315497_wp_
+    b = f * a - b +   0.000192469359688_wp_
+    a = f * b - a +   0.003416331766012_wp_
+    b = f * a - b +   0.047718748798174_wp_
+    a = f * b - a +   0.509493365439983_wp_
+    b = f * a - b +   4.011673760179349_wp_
+    a = f * b - a +  22.274819242462231_wp_
+    b = f * a - b +  82.489032744024100_wp_
+    a = f * b - a + 190.494320172742844_wp_
+    a = f * a - b + 255.466879624362167_wp_
+    besik=0.5_wp_*(a-b)
+    if(l .and. e) besik=exp(-v)*besik
+    if(l) return
+!
+    a =           +   0.000000000000003_wp_
+    b = f * a     +   0.000000000000159_wp_
+    a = f * b - a +   0.000000000007658_wp_
+    b = f * a - b +   0.000000000318588_wp_
+    a = f * b - a +   0.000000011281211_wp_
+    b = f * a - b +   0.000000335195256_wp_
+    a = f * b - a +   0.000008216025940_wp_
+    b = f * a - b +   0.000162708379043_wp_
+    a = f * b - a +   0.002536308188086_wp_
+    b = f * a - b +   0.030080722420512_wp_
+    a = f * b - a +   0.259084432434900_wp_
+    b = f * a - b +   1.511535676029228_wp_
+    a = f * b - a +   5.283632866873920_wp_
+    b = f * a - b +   8.005368868700334_wp_
+    a = f * b - a -   4.563433586448395_wp_
+    a = f * a - b -  21.057660177402440_wp_
+    besik=-log(0.125_wp_*x)*besik+0.5_wp_*(a-b)
+    if(e) besik=exp(x)*besik
+    return
+!
+      4 f=32.0_wp_/v-2.0_wp_
+    b =           - 0.000000000000001_wp_
+    a = f * b     - 0.000000000000001_wp_
+    b = f * a - b + 0.000000000000004_wp_
+    a = f * b - a + 0.000000000000010_wp_
+    b = f * a - b - 0.000000000000024_wp_
+    a = f * b - a - 0.000000000000104_wp_
+    b = f * a - b + 0.000000000000039_wp_
+    a = f * b - a + 0.000000000000966_wp_
+    b = f * a - b + 0.000000000001800_wp_
+    a = f * b - a - 0.000000000004497_wp_
+    b = f * a - b - 0.000000000033127_wp_
+    a = f * b - a - 0.000000000078957_wp_
+    b = f * a - b + 0.000000000029802_wp_
+    a = f * b - a + 0.000000001238425_wp_
+    b = f * a - b + 0.000000008513091_wp_
+    a = f * b - a + 0.000000056816966_wp_
+    b = f * a - b + 0.000000513587727_wp_
+    a = f * b - a + 0.000007247591100_wp_
+    b = f * a - b + 0.000172700630778_wp_
+    a = f * b - a + 0.008445122624921_wp_
+    a = f * a - b + 2.016558410917480_wp_
+    besik=0.199471140200717_wp_*(a-b)/sqrt(v)
+    if(e) return
+    besik=exp(v)*besik
+    return
+!
+    entry ebesi1(x)
+!
+    e=.true.
+    go to 2
+!
+    entry besi1(x)
+!
+    e=.false.
+      2 l=.true.
+    v=abs(x)
+    if(v >= 8.0_wp_) go to 3
+      7 f=0.0625_wp_*x**2-2.0_wp_
+    a =           +   0.000000000000001_wp_
+    b = f * a     +   0.000000000000031_wp_
+    a = f * b - a +   0.000000000001679_wp_
+    b = f * a - b +   0.000000000079291_wp_
+    a = f * b - a +   0.000000003227617_wp_
+    b = f * a - b +   0.000000111946285_wp_
+    a = f * b - a +   0.000003264138122_wp_
+    b = f * a - b +   0.000078756785754_wp_
+    a = f * b - a +   0.001543019015627_wp_
+    b = f * a - b +   0.023993079147841_wp_
+    a = f * b - a +   0.287855511804672_wp_
+    b = f * a - b +   2.571459906347755_wp_
+    a = f * b - a +  16.334550552522066_wp_
+    b = f * a - b +  69.395917633734448_wp_
+    a = f * b - a + 181.312616040570265_wp_
+    a = f * a - b + 259.890237806477292_wp_
+    besik=0.0625_wp_*(a-b)*x
+    if(l .and. e) besik=exp(-v)*besik
+    if(l) return
+!
+    a =           +   0.000000000000001_wp_
+    b = f * a     +   0.000000000000042_wp_
+    a = f * b - a +   0.000000000002163_wp_
+    b = f * a - b +   0.000000000096660_wp_
+    a = f * b - a +   0.000000003696783_wp_
+    b = f * a - b +   0.000000119367971_wp_
+    a = f * b - a +   0.000003202510692_wp_
+    b = f * a - b +   0.000070010627855_wp_
+    a = f * b - a +   0.001217056994516_wp_
+    b = f * a - b +   0.016300049289816_wp_
+    a = f * b - a +   0.161074301656148_wp_
+    b = f * a - b +   1.101461993004852_wp_
+    a = f * b - a +   4.666387026862842_wp_
+    b = f * a - b +   9.361617831395389_wp_
+    a = f * b - a -   1.839239224286199_wp_
+    a = f * a - b -  26.688095480862668_wp_
+    besik=log(0.125_wp_*x)*besik+1.0_wp_/x-0.0625_wp_*(a-b)*x
+    if(e) besik=exp(x)*besik
+    return
+!
+      3 f=32.0_wp_/v-2.0_wp_
+    b =           + 0.000000000000001_wp_
+    a = f * b     + 0.000000000000001_wp_
+    b = f * a - b - 0.000000000000005_wp_
+    a = f * b - a - 0.000000000000010_wp_
+    b = f * a - b + 0.000000000000026_wp_
+    a = f * b - a + 0.000000000000107_wp_
+    b = f * a - b - 0.000000000000053_wp_
+    a = f * b - a - 0.000000000001024_wp_
+    b = f * a - b - 0.000000000001804_wp_
+    a = f * b - a + 0.000000000005103_wp_
+    b = f * a - b + 0.000000000035408_wp_
+    a = f * b - a + 0.000000000081531_wp_
+    b = f * a - b - 0.000000000047563_wp_
+    a = f * b - a - 0.000000001401141_wp_
+    b = f * a - b - 0.000000009613873_wp_
+    a = f * b - a - 0.000000065961142_wp_
+    b = f * a - b - 0.000000629724239_wp_
+    a = f * b - a - 0.000009732146728_wp_
+    b = f * a - b - 0.000277205360764_wp_
+    a = f * b - a - 0.024467442963276_wp_
+    a = f * a - b + 1.951601204652572_wp_
+    besik=0.199471140200717_wp_*(a-b)/sqrt(v)
+    if(x < 0.0_wp_) besik=-besik
+    if(e) return
+    besik=exp(v)*besik
+    return
+!
+    entry ebesk0 (x)
+!
+    e=.true.
+    go to 11
+!
+    entry besk0(x)
+!
+    e=.false.
+     11 if(x <= 0.0_wp_) go to 9
+    l=.false.
+    v=x
+    if(x < 5.0_wp_) go to 8
+    f=20.0_wp_/x-2.0_wp_
+    a =           - 0.000000000000002_wp_
+    b = f * a     + 0.000000000000011_wp_
+    a = f * b - a - 0.000000000000079_wp_
+    b = f * a - b + 0.000000000000581_wp_
+    a = f * b - a - 0.000000000004580_wp_
+    b = f * a - b + 0.000000000039044_wp_
+    a = f * b - a - 0.000000000364547_wp_
+    b = f * a - b + 0.000000003792996_wp_
+    a = f * b - a - 0.000000045047338_wp_
+    b = f * a - b + 0.000000632575109_wp_
+    a = f * b - a - 0.000011106685197_wp_
+    b = f * a - b + 0.000269532612763_wp_
+    a = f * b - a - 0.011310504646928_wp_
+    a = f * a - b + 1.976816348461652_wp_
+    besik=0.626657068657750_wp_*(a-b)/sqrt(x)
+    if(e) return
+    z=besik
+    besik=0.0_wp_
+    if(x < 180.0_wp_) besik=exp(-x)*z
+    return
+!
+    entry ebesk1(x)
+!
+    e=.true.
+    go to 12
+!
+    entry besk1(x)
+!
+    e=.false.
+     12 if(x <= 0.0_wp_) go to 9
+    l=.false.
+    v=x
+    if(x < 5.0_wp_) go to 7
+    f=20.0_wp_/x-2.0_wp_
+    a =           + 0.000000000000002_wp_
+    b = f * a     - 0.000000000000013_wp_
+    a = f * b - a + 0.000000000000089_wp_
+    b = f * a - b - 0.000000000000663_wp_
+    a = f * b - a + 0.000000000005288_wp_
+    b = f * a - b - 0.000000000045757_wp_
+    a = f * b - a + 0.000000000435417_wp_
+    b = f * a - b - 0.000000004645555_wp_
+    a = f * b - a + 0.000000057132218_wp_
+    b = f * a - b - 0.000000845172048_wp_
+    a = f * b - a + 0.000016185063810_wp_
+    b = f * a - b - 0.000468475028167_wp_
+    a = f * b - a + 0.035465291243331_wp_
+    a = f * a - b + 2.071901717544716_wp_
+    besik=0.626657068657750_wp_*(a-b)/sqrt(x)
+    if(e) return
+    z=besik
+    besik=0.0_wp_
+    if(x < 180.0_wp_) besik=exp(-x)*z
+    return
+      9 besik=0.0_wp_
+    write(nout,"(1x,'besik ... non-positive argument x = ',e15.4)") x
+  end function besik
+!
+! routines for conical function: end
+!
+end module conical

src/const_and_precisions.f90

@@ -1,17 +1,21 @@
 !########################################################################!
 
  MODULE const_and_precisions
- use itm_types, only : wp_ => r8
- use itm_constants, only : pi => itm_pi, e_ => itm_qe, me_ => itm_me, c_ => itm_c
 !########################################################################!
  IMPLICIT NONE
  PUBLIC
 !------------------------------------------------------------------------
 ! common precisions
 !------------------------------------------------------------------------
-! INTEGER, PARAMETER :: sp_   = 4         ! single precision
-! INTEGER, PARAMETER :: dp_   = 8         ! double precision
-! INTEGER, PARAMETER :: wp_   = dp_       ! work-precision
+!  INTEGER, PARAMETER :: i1 = SELECTED_INT_KIND (2)        ! Integer*1
+!  INTEGER, PARAMETER :: i2 = SELECTED_INT_KIND (4)        ! Integer*2
+  INTEGER, PARAMETER :: i4 = SELECTED_INT_KIND (9)        ! Integer*4
+  INTEGER, PARAMETER :: i8 = SELECTED_INT_KIND (18)       ! Integer*8
+  INTEGER, PARAMETER :: r4 = SELECTED_REAL_KIND (6, 37)   ! Real*4
+  INTEGER, PARAMETER :: r8 = SELECTED_REAL_KIND (15, 300) ! Real*8
+! INTEGER, PARAMETER :: sp_   = r4        ! single precision
+! INTEGER, PARAMETER :: dp_   = r8        ! double precision
+  INTEGER, PARAMETER :: wp_   = r8        ! work-precision
 ! INTEGER, PARAMETER :: odep_ = dp_       ! ODE-solver precision
 ! INTEGER, PARAMETER :: xp_   = wp_       ! for ext. modules if necessary
 !------------------------------------------------------------------------
@@ -26,31 +30,36 @@
 !!========================================================================
 ! Arithmetic constants
 !========================================================================
+  integer,   parameter :: izero   = 0
   REAL(wp_), PARAMETER :: zero    = 0.0_wp_
-  REAL(wp_), PARAMETER :: unit    = 1.0_wp_
-!  REAL(wp_), PARAMETER :: pi      = 3.141592653589793_wp_
-! REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
+  REAL(wp_), PARAMETER :: half    = 0.5_wp_
+  REAL(wp_), PARAMETER :: one     = 1.0_wp_
+  REAL(wp_), PARAMETER :: two     = 2.0_wp_
+  real(wp_), parameter :: pi      = 3.141592653589793_wp_  ! 3.141592653589793238462643383280
+  real(wp_), parameter :: pihalf  = 1.57079632679489661923_wp_
+  REAL(wp_), PARAMETER :: sqrt_pi = 1.772453850905516_wp_
 ! REAL(wp_), PARAMETER :: sqrt_2  = 1.414213562373095_wp_
-! REAL(wp_), PARAMETER :: rad     = pi/180.0_wp_
+  REAL(wp_), PARAMETER :: degree  = pi/180.0_wp_
+  REAL(wp_), PARAMETER :: emn1    = 0.367879441171442_wp_  ! exp(-1)
 !---
-! REAL(wp_), PARAMETER :: ex(1:3) = (/unit,zero,zero/)
-! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,unit,zero/)
-! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,unit/)
+! REAL(wp_), PARAMETER :: ex(1:3) = (/one ,zero,zero/)
+! REAL(wp_), PARAMETER :: ey(1:3) = (/zero,one ,zero/)
+! REAL(wp_), PARAMETER :: ez(1:3) = (/zero,zero,one /)
 !---
-! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/unit,zero,zero, &
-!                                               zero,unit,zero, &
-!                                               zero,zero,unit/),(/3,3/))
-! COMPLEX(wp_), PARAMETER :: im    = (0.0_wp_,1.0_wp_)
-! COMPLEX(wp_), PARAMETER :: czero = (0.0_wp_,0.0_wp_)
-! COMPLEX(wp_), PARAMETER :: cunit = (1.0_wp_,0.0_wp_)
+! REAL(wp_), PARAMETER :: kron(3,3) = reshape((/one ,zero,zero, &
+!                                               zero,one ,zero, &
+!                                               zero,zero,one /),(/3,3/))
+  COMPLEX(wp_), PARAMETER :: im    = (0.0_wp_,1.0_wp_)
+  COMPLEX(wp_), PARAMETER :: czero = (0.0_wp_,0.0_wp_)
+  COMPLEX(wp_), PARAMETER :: cunit = (1.0_wp_,0.0_wp_)
 ! COMPLEX(wp_), PARAMETER :: ctwo  = (2.0_wp_,0.0_wp_)
 !========================================================================
 ! Computer constants
 !========================================================================
-  REAL(wp_), PARAMETER :: comp_eps      = EPSILON(unit)
-!  REAL(wp_), PARAMETER :: comp_eps2     = comp_eps**2
-! REAL(wp_), PARAMETER :: comp_tiny     = TINY(unit)
-! REAL(wp_), PARAMETER :: comp_huge     = HUGE(unit)
+  REAL(wp_), PARAMETER :: comp_eps      = EPSILON(one)
+! REAL(wp_), PARAMETER :: comp_eps2     = comp_eps**2
+  REAL(wp_), PARAMETER :: comp_tiny     = TINY(one)
+  REAL(wp_), PARAMETER :: comp_huge     = HUGE(one)
 ! REAL(wp_), PARAMETER :: comp_tinylog  =-200 ! LOG10(comp_tiny)
 ! REAL(wp_), PARAMETER :: comp_hugelog  =+200 ! LOG10(comp_huge)
 ! REAL(wp_), PARAMETER :: comp_tiny1    = 1d+50*comp_tiny
@@ -60,26 +69,42 @@
 !------------------------------------------------------------------------
 ! Conventional constants
 !------------------------------------------------------------------------
+  INTEGER,  PARAMETER :: int_invalid = -999999999
+  REAL(R8), PARAMETER :: r8_invalid = -9.0e40_r8
 ! REAL(wp_), PARAMETER :: output_tiny = 1.0d-66
 ! REAL(wp_), PARAMETER :: output_huge = 1.0d+66
 !========================================================================
 ! Physical constants (SI)
 !========================================================================
-!  REAL(wp_), PARAMETER :: e_     = 1.602176487d-19    ! [C]
-!  REAL(wp_), PARAMETER :: me_    = 9.10938215d-31    ! [kg]
-! REAL(wp_), PARAMETER :: mp_    = 1.672621637d-27    ! [kg]
-! REAL(wp_), PARAMETER :: rmpe_  = mp_/me_
-!  REAL(wp_), PARAMETER :: c_     = 2.99792458d+08    ! [m/s]
-! REAL(wp_), PARAMETER :: eps0_  = 8.854187817d-12    ! [F/m]
+  real (wp_), parameter :: e_ = 1.602176487e-19_wp_         ! elementary charge, C
+  real (wp_), parameter :: me_ = 9.10938215e-31_wp_         ! electron mass, kg
+! real (wp_), parameter :: mp_ = 1.672621637e-27_wp_        ! proton mass, kg
+! real (wp_), parameter :: md_ = 3.34358320e-27_wp_         ! deuteron mass, kg
+! real (wp_), parameter :: mt_ = 5.00735588e-27_wp_         ! triton mass, kg
+! real (wp_), parameter :: ma_ = 6.64465620e-27_wp_         ! alpha mass, kg
+! real (wp_), parameter :: amu_ = 1.660538782e-27_wp_       ! amu, kg
+! REAL (wp_), PARAMETER :: rmpe_  = mp_/me_                 ! proton-electron mass ratio
+  real (wp_), parameter :: c_ = 2.99792458e8_wp_            ! speed of light, m/s
+  real (wp_), parameter :: mu0_ = 4.0e-7_wp_ * pi           ! magnetic permeability of vacuum
+  real (wp_), parameter :: eps0_ = 1.0_wp_ / (mu0_ * c_**2) ! dielectric constant of vacuum, F/m
+! real (wp_), parameter :: avogr = 6.02214179e23_wp_
+! real (wp_), parameter :: KBolt = 1.3806504e-23_wp_
+!========================================================================
+! Physical constants (cgs)
+!========================================================================
+  real (wp_), parameter :: ccgs_ = c_*1.e2_wp_              ! speed of light, cm/s
+  real (wp_), parameter :: mecgs_ = me_*1.e3_wp_            ! electron mass, g
+  real (wp_), parameter :: ecgs_ = e_*c_*10._wp_            ! elementary charge, statcoul
 !------------------------------------------------------------------------
 ! Useful definitions
 !------------------------------------------------------------------------
-  REAL(wp_), PARAMETER :: keV_   = 1000*e_         ! [J]
+  REAL(wp_), PARAMETER :: keV_   = 1.e3_wp_*e_     ! [J]
   REAL(wp_), PARAMETER :: mc2_SI = me_*c_**2       ! [J]
   REAL(wp_), PARAMETER :: mc2_   = mc2_SI/keV_     ! [keV]
+  REAL(wp_), PARAMETER :: mu0inv = 1._wp_/mu0_     !
 ! REAL(wp_), PARAMETER :: mc_    = me_*c_          ! [kg*m/s]
 ! ! f_ce = fce1_*B (B in Tesla):                   !
-! REAL(wp_), PARAMETER :: wce1_  = e_/me_          ! [rad/s]
+  REAL(wp_), PARAMETER :: wce1_  = e_/me_          ! [rad/s]
 ! REAL(wp_), PARAMETER :: fce1_  = wce1_/(2*pi)    ! [1/s]
 ! ! f_pl = fpe1_*sqrt(Ne) (Ne in 1/m**3):          !
 ! REAL(wp_), PARAMETER :: wpe1_  = 56.4049201      ! [rad/s]
@@ -100,6 +125,33 @@
 ! REAL(wp_), PARAMETER :: Npar_min = 1.0d-3
 !########################################################################!
 
+  interface is_valid
+     module procedure is_valid_int4, is_valid_int8, is_valid_real8
+  end interface
+
+contains
+
+  logical function is_valid_int4(in_int)
+    implicit none
+    integer(i4), intent(in) :: in_int
+    is_valid_int4 = in_int /= int_invalid
+    return
+  end function is_valid_int4
+
+  logical function is_valid_int8(in_int)
+    implicit none
+    integer(i8), intent(in) :: in_int
+    is_valid_int8 = in_int /= int_invalid
+    return
+  end function is_valid_int8
+
+  logical function is_valid_real8(in_real)
+    implicit none
+    real(r8), intent(in) :: in_real
+    is_valid_real8 = abs(in_real - r8_invalid) > abs(r8_invalid) * 1.0e-15_r8
+    return
+  end function is_valid_real8
+
  END MODULE const_and_precisions
 
 !########################################################################!

src/coreprofiles.f90

@@ -0,0 +1,319 @@
+module coreprofiles
+  use const_and_precisions, only : wp_,zero,one
+  implicit none
+
+  INTEGER, SAVE :: npp,nsfd
+  REAL(wp_), SAVE :: psdbnd,psnpp,denpp,ddenpp,d2denpp
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfn,cfn,psrad
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: ct,cz
+  REAL(wp_), SAVE :: dens0,aln1,aln2,te0,dte0,alt1,alt2,zeffan
+
+contains
+
+  subroutine density(psin,dens,ddens)
+    use graydata_flags, only : iprof
+!    use graydata_anequil, only : dens0,aln1,aln2
+    use dierckx, only : splev,splder
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_), intent(out) :: dens,ddens
+! local variables
+    integer, parameter :: nn=3, nn1=nn+1, nn2=nn+2
+    integer :: ier,nu
+    real(wp_) :: profd,dprofd,dpsib,tt,fp,dfp,fh,dfh
+    real(wp_), dimension(1) :: xxs,ffs
+    real(wp_), dimension(npp+4) :: wrkfd
+
+!
+!     computation of density [10^19 m^-3] and derivative wrt psi
+!
+    dens=zero
+    ddens=zero
+    if((psin >= psdbnd).or.(psin < zero)) return
+!
+    if(iprof == 0) then
+      if(psin > one) return
+      profd=(one-psin**aln1)**aln2
+      dens=dens0*profd
+      dprofd=-aln1*aln2*psin**(aln1-one)           &
+              *(one-psin**aln1)**(aln2-one)
+      ddens=dens0*dprofd
+    else
+      if(psin > psnpp) then
+
+!     smooth interpolation for psnpp < psi < psdbnd
+!     dens = fp * fh
+!     fp: parabola matched at psi=psnpp with given profile density
+!     fh=(1-t)^3(1+3t+6t^2) is a smoothing function: 
+!     fh(0)=1, fh(1)=0 and zero first and second deriv at t=0,1 
+!
+        dpsib=psin-psnpp
+        fp=denpp+dpsib*ddenpp+0.5_wp_*dpsib**2*d2denpp
+        dfp=ddenpp+dpsib*d2denpp
+        tt=dpsib/(psdbnd-psnpp)
+        fh=(one-tt)**3*(one+3.0_wp_*tt+6.0_wp_*tt**2)
+        dfh=-30.0_wp_*(one-tt)**2*tt**2/(psdbnd-psnpp)
+        dens=fp*fh
+        ddens=dfp*fh+fp*dfh
+      else
+        xxs(1)=psin
+        ier=0
+        call splev(tfn,nsfd,cfn,3,xxs,ffs,1,ier)
+        dens=ffs(1)
+        nu=1
+        ier=0
+        call splder(tfn,nsfd,cfn,3,nu,xxs,ffs,1,wrkfd,ier)
+        ddens=ffs(1)
+        if(ier > 0) print*,ier
+        if(abs(dens) < 1.0e-10_wp_) dens=zero
+      end if
+!      if(dens < zero) print*,'    DENSITY NEGATIVE',dens
+      if(dens < zero) then
+        dens=zero
+        ddens=zero
+      end if  
+    end if
+  end subroutine density
+
+  function temp(psin)
+    use const_and_precisions, only : wp_,zero,one
+    use graydata_flags, only : iprof
+!    use graydata_anequil, only : te0,dte0,alt1,alt2
+    use utils, only : locate
+    use simplespline, only :spli
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_) :: temp
+! local variables
+    integer :: k
+    real(wp_) :: proft,dps
+
+    temp=zero
+    if((psin >= one).or.(psin < zero)) return
+    if(iprof == 0) then
+      proft=(1.0_wp_-psin**alt1)**alt2
+      temp=(te0-dte0)*proft+dte0
+    else
+      call locate(psrad,npp,psin,k)
+      k=max(1,min(k,npp-1))
+      dps=psin-psrad(k)
+      temp=spli(ct,npp,k,dps)
+    endif
+  end function temp
+
+  function fzeff(psin)
+    use const_and_precisions, only : wp_,zero,one
+    use graydata_flags, only : iprof
+!    use graydata_anequil, only : zeffan
+    use utils, only : locate
+    use simplespline, only :spli
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_) :: fzeff
+! local variables
+    integer :: k
+    real(wp_) :: dps
+
+    fzeff=one
+    if((psin >= one).or.(psin < zero)) return
+    if(iprof == 0) then
+      fzeff=zeffan
+    else
+      call locate(psrad,npp,psin,k)
+      k=max(1,min(k,npp-1))
+      dps=psin-psrad(k)
+      fzeff=spli(cz,npp,k,dps)
+    endif
+  end function fzeff
+  
+  subroutine read_profiles(filenm,psin,te,ne,zeff,unit)
+    use utils, only : get_free_unit
+    implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), dimension(:), allocatable, intent(out) :: psin,te,ne,zeff
+    integer, optional, intent(in) :: unit
+! local variables
+    integer :: u, i, n
+    
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+    open(file=trim(filenm),status='old',unit=u)
+    read(u,*) n
+    if(allocated(psin)) deallocate(psin)
+    if(allocated(te)) deallocate(te)
+    if(allocated(ne)) deallocate(ne)
+    if(allocated(zeff)) deallocate(zeff)
+    allocate(psin(n),te(n),ne(n),zeff(n))
+    do i=1,n
+      read(u,*) psin(i),te(i),ne(i),zeff(i)
+    end do
+    psin(1)=max(psin(1),zero)
+    close(u)
+  end subroutine read_profiles
+
+  subroutine read_profiles_an(filenm,te,ne,zeff,unit)
+    use utils, only : get_free_unit
+     implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), dimension(:), allocatable, intent(out) :: te,ne,zeff
+    integer, optional, intent(in) :: unit
+! local variables
+    integer :: u
+        
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+
+    if(allocated(te)) deallocate(te)
+    if(allocated(ne)) deallocate(ne)
+    if(allocated(zeff)) deallocate(zeff)
+    allocate(te(4),ne(3),zeff(1))
+
+    open(file=trim(filenm),status='old',unit=u)
+    read(u,*) ne(1:3) ! dens0,aln1,aln2
+    read(u,*) te(1:4) ! te0,dte0,alt1,alt2
+    read(u,*) zeff(1) ! zeffan
+    close(u)
+  end subroutine read_profiles_an
+
+  subroutine tene_scal(te,ne,tfact,nfact,bfact,iscal)
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(inout) :: te,ne
+    real(wp_), intent(in) :: tfact,nfact,bfact
+    integer, intent(in) :: iscal
+! local variables
+    real(wp_) :: aat,aan,ffact
+
+    if (iscal==0) then
+      aat=2.0_wp_/3.0_wp_
+      aan=4.0_wp_/3.0_wp_
+    else
+      aat=1.0_wp_
+      aan=1.0_wp_
+    end if
+    if(iscal==2) then
+      ffact=1.0_wp_
+    else
+      ffact=bfact
+    end if
+    te(:)=te(:)*ffact**aat*tfact
+    ne(:)=ne(:)*ffact**aan*nfact
+  end subroutine tene_scal
+
+  subroutine set_prfspl(psin,te,ne,zeff,ssplne)
+    use simplespline, only : difcs
+    use dierckx, only : curfit, splev, splder
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(in) :: psin,te,ne,zeff
+    real(wp_), intent(in) :: ssplne
+! local variables
+    integer, parameter :: iopt=0, kspl=3
+    integer :: n, npest, lwrkf, ier
+    real(wp_) :: xb, xe, fp, xnv, ynv
+    real(wp_), dimension(:), allocatable :: wf, wrkf
+    integer, dimension(:), allocatable :: iwrkf
+    real(wp_), dimension(1) :: dedge,ddedge,d2dedge
+
+    n=size(psin)
+    npest=n+4
+    lwrkf=n*4+npest*16
+    allocate(wrkf(lwrkf),iwrkf(npest),wf(n))
+
+!   if necessary, reallocate spline arrays
+    if(.not.allocated(psrad)) then
+      allocate(psrad(n),ct(n,4),cz(n,4))
+    else
+      if(size(psrad)<n) then
+        deallocate(psrad,ct,cz)
+        allocate(psrad(n),ct(n,4),cz(n,4))
+      end if
+    end if
+    if(.not.allocated(cfn)) then
+      allocate(tfn(npest),cfn(npest))
+    else
+      if(size(cfn)<npest) then
+        deallocate(tfn,cfn)
+        allocate(tfn(npest),cfn(npest))
+      end if
+    end if
+
+!   spline approximation of  temperature and Zeff
+    call difcs(psin,te,  n,iopt,ct,ier)
+    call difcs(psin,zeff,n,iopt,cz,ier)
+    psrad=psin
+    npp=n
+
+!   spline approximation of density
+    xb=zero
+    xe=psin(n)
+    wf(:)=one
+    call curfit(iopt,n,psin,ne,wf,xb,xe,kspl,ssplne,npest, &
+       nsfd,tfn,cfn,fp,wrkf,lwrkf,iwrkf,ier)
+
+!   compute polinomial extrapolation matching the spline boundary up to the
+!   2nd order derivative, extending the profile up to psi=psdbnd where
+!   ne=ne'=ne''=0
+!   spline value and derivatives at the edge
+    call splev(tfn,nsfd,cfn,3,psin(n:n),dedge(1:1),1,ier)
+    call splder(tfn,nsfd,cfn,3,1,psin(n:n),ddedge(1:1), 1,wrkf(1:nsfd),ier)
+    call splder(tfn,nsfd,cfn,3,2,psin(n:n),d2dedge(1:1),1,wrkf(1:nsfd),ier)
+!   determination of the boundary
+    psnpp=psin(n)
+    denpp=dedge(1)
+    ddenpp=ddedge(1)
+    d2denpp=d2dedge(1)
+    psdbnd=psnpp
+    if(ddenpp < zero) then
+      xnv=psnpp-ddenpp/d2denpp
+      ynv=denpp-0.5_wp_*ddenpp**2/d2denpp
+      if((d2denpp > zero).and.(ynv >= zero))  then
+        psdbnd=xnv
+      else
+        psdbnd=xnv+sqrt((ddenpp/d2denpp)**2-2.0_wp_*denpp/d2denpp)
+      end if  
+      print*,psnpp,denpp,ddenpp,d2denpp,xnv,ynv
+      print*,psdbnd
+    end if
+
+    deallocate(iwrkf,wrkf,wf)
+  end subroutine set_prfspl
+
+  subroutine unset_prfspl
+    implicit none
+
+    if(allocated(psrad)) deallocate(psrad)
+    if(allocated(ct)) deallocate(ct)
+    if(allocated(cz)) deallocate(cz)
+    if(allocated(tfn)) deallocate(tfn)
+    if(allocated(cfn)) deallocate(cfn)
+  end subroutine unset_prfspl
+
+  subroutine set_prfan(te,ne,zeff)
+    implicit none
+    REAL(wp_), dimension(:), intent(in) :: te,ne,zeff
+
+    te0=te(1)
+    dte0=te(2)
+    alt1=te(3)
+    alt2=te(4)
+    dens0=ne(1)
+    aln1=ne(2)
+    aln2=ne(3)
+    zeffan=zeff(1)
+
+    psdbnd=1.0_wp_
+  end subroutine set_prfan
+  
+end module coreprofiles

src/eccd.f90

@@ -0,0 +1,884 @@
+module eccd
+  use const_and_precisions, only : wp_
+  implicit none
+  real(wp_), parameter, private :: cst2min=1.0e-6_wp_ ! min width of trap. cone
+  integer, parameter, private :: nfpp=13, & ! number of extra parameters passed
+            nfpp1=nfpp+ 1, nfpp2=nfpp+ 2, & ! to the integrand function fpp
+            nfpp3=nfpp+ 3, nfpp4=nfpp+ 4, &
+            nfpp5=nfpp+ 5
+!########################################################################
+! the following parameters are used by N.M. subroutines:
+! The module contains few subroutines which are requested to calculate 
+! the current drive value by adjoint approach
+!########################################################################
+   CHARACTER, PRIVATE, PARAMETER :: adj_appr(1:6) = & ! adj. approach switcher
+                    (/ 'l',  &     ! (1)='l': collisionless limit
+                                   ! (1)='c': collisional (classical) limit,
+                                   !          w/o trap. part.
+                       'm',  &     ! (2)='m': momentum conservation
+                                   ! (2)='h': high-speed limit
+!--- 
+                       'l',  &     ! DO NOT CHANGE!
+                       'r',  &     ! DO NOT CHANGE!
+                       'v',  &     ! DO NOT CHANGE!
+                       'i' /)      ! DO NOT CHANGE!
+!-------
+   REAL(wp_), PRIVATE :: r2,q2,gp1 ! coefficients for HSL integrand function
+!-------
+   REAL(wp_), PRIVATE, PARAMETER :: delta = 1e-4   ! border for recalculation
+!------- for N.M. subroutines (variational principle) -------
+   REAL(wp_), PRIVATE :: sfd(1:4)  ! polyn. exp. of the "Spitzer"-function
+   INTEGER,   PRIVATE, PARAMETER :: nre = 2        ! order of rel. correct.
+   REAL(wp_), PRIVATE, PARAMETER :: vp_mee(0:4,0:4,0:2) = &
+            RESHAPE((/0.0,      0.0,      0.0,      0.0,      0.0,      &
+                      0.0,      0.184875, 0.484304, 1.06069,  2.26175,  &
+                      0.0,      0.484304, 1.41421,  3.38514,  7.77817,  &
+                      0.0,      1.06069,  3.38514,  8.73232,  21.4005,  &
+                      0.0,      2.26175,  7.77817,  21.4005,  55.5079,  &
+                      !                                                 &
+                      0.0,     -1.33059,-2.57431, -5.07771, -10.3884,   &
+                     -0.846284,-1.46337, -1.4941, -0.799288, 2.57505,   &
+                     -1.1601,  -1.4941,   2.25114,  14.159,   50.0534,  &
+                     -1.69257, -0.799288, 14.159,   61.4168,  204.389,  &
+                     -2.61022,  2.57505,  50.0534,  204.389,  683.756,  &
+                      !                                                 &
+                      0.0,      2.62498,  0.985392,-5.57449, -27.683,   &
+                      0.0,      3.45785,  5.10096,  9.34463,  22.9831,  &
+                     -0.652555, 5.10096,  20.5135,  75.8022,  268.944,  &
+                     -2.11571,  9.34463,  75.8022,  330.42,   1248.69,  &
+                     -5.38358,  22.9831,  268.944,  1248.69,  4876.48/),&
+                    (/5,5,3/))
+   REAL(wp_), PRIVATE, PARAMETER :: vp_mei(0:4,0:4,0:2) = &
+            RESHAPE((/0.0,     0.886227, 1.0,      1.32934,  2.0,       &
+                      0.886227,1.0,      1.32934,  2.0,      3.32335,   &
+                      1.0,     1.32934,  2.0,      3.32335,  6.0,       &
+                      1.32934, 2.0,      3.32335,  6.0,      11.6317,   &
+                      2.0,     3.32335,  6.0,      11.6317,  24.0,      &
+                      !                                                 &
+                      0.0,     0.332335, 1.0,      2.49251,  6.0,       &
+                      1.66168, 1.0,      2.49251,  6.0,      14.5397,   &
+                      3.0,     2.49251,  6.0,      14.5397,  36.0,      &
+                      5.81586, 6.0,      14.5397,  36.0,     91.5999,   &
+                      12.0,    14.5397,  36.0,     91.5999,  240.0,     &
+                      !                                                 &
+                      0.0,    -0.103855, 0.0,      1.09047,  6.0,       &
+                      0.726983,0.0,      1.09047,  6.0,      24.5357,   &
+                      3.0,     1.09047,  6.0,      24.5357,  90.0,      &
+                      9.81427, 6.0,      24.5357,  90.0,     314.875,   &
+                      30.0,    24.5357,  90.0,     314.875,  1080.0 /), &
+                    (/5,5,3/))
+   REAL(wp_), PRIVATE, PARAMETER :: vp_oee(0:4,0:4,0:2) = &
+            RESHAPE((/0.0,     0.56419,  0.707107, 1.0073,   1.59099,   &
+                      0.56419, 0.707107, 1.0073,   1.59099,  2.73981,   &
+                      0.707107,1.0073,   1.59099,  2.73981,  5.08233,   &
+                      1.0073,  1.59099,  2.73981,  5.08233,  10.0627,   &
+                      1.59099, 2.73981,  5.08233,  10.0627,  21.1138,   &
+                      !                                                 &
+                      0.0,     1.16832,  1.90035,  3.5758,   7.41357,   &
+                      2.17562, 1.90035,  3.5758,   7.41357,  16.4891,   &
+                      3.49134, 3.5758,   7.41357,  16.4891,  38.7611,   &
+                      6.31562, 7.41357,  16.4891,  38.7611,  95.4472,   &
+                      12.4959, 16.4891,  38.7611,  95.4472,  244.803,   &
+                      !                                                 &
+                      0.0,     2.65931,  4.64177,  9.6032,   22.6941,   &
+                      4.8652,  4.64177,  9.6032,   22.6941,  59.1437,   &
+                      9.51418, 9.6032,   22.6941,  59.1437,  165.282,   &
+                      21.061,  22.6941,  59.1437,  165.282,  485.785,   &
+                      50.8982, 59.1437,  165.282,  485.785,  1483.22/), &
+                    (/5,5,3/))
+   REAL(wp_), PRIVATE, PARAMETER :: vp_g(0:4,0:2) = &
+            RESHAPE((/1.32934, 2.0,      3.32335,  6.0,      11.6317,   &
+                      2.49251, 0.0,      2.90793,  12.0,     39.2571,   &
+                      1.09047, 6.0,      11.45,    30.0,     98.9606/), &
+                    (/5,3/))
+!########################################################################
+
+  interface setcdcoeff
+    module procedure setcdcoeff_notrap,setcdcoeff_cohen,setcdcoeff_ncl
+  end interface setcdcoeff
+
+contains
+
+  subroutine setcdcoeff_notrap(zeff,cst2,eccdpar)
+    implicit none
+    real(wp_), intent(in) :: zeff
+    real(wp_), intent(out) :: cst2
+    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+
+    cst2=0.0_wp_
+    allocate(eccdpar(1))
+    eccdpar(1)=zeff
+  end subroutine setcdcoeff_notrap
+
+  subroutine setcdcoeff_cohen(zeff,rbn,rbx,cst2,eccdpar)
+!     cohen model      
+!     rbn=B/B_min
+!     rbx=B/B_max
+!     cst2=1.0_wp_-B/B_max
+!     alams=sqrt(1-B_min/B_max)
+!     Zeff < 31 !!!
+!     fp0s= P_a (alams)
+    use conical, only : fconic
+    implicit none
+    real(wp_), intent(in) :: zeff,rbn,rbx
+    real(wp_), intent(out) :: cst2
+    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+    real(wp_) :: alams,pa,fp0s
+
+    cst2=1.0_wp_-rbx
+    if(cst2<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)
+    allocate(eccdpar(5))
+    eccdpar(1)=zeff
+    eccdpar(2)=rbn
+    eccdpar(3)=alams
+    eccdpar(4)=pa
+    eccdpar(5)=fp0s
+  end subroutine setcdcoeff_cohen
+
+  subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cst2,eccdpar)
+    use magsurf_data, only : ch,tjp,tlm,njpt,nlmt
+    use dierckx, only : profil
+    implicit none
+    integer, parameter :: ksp=3
+    real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
+    real(wp_), intent(out) :: cst2
+    real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
+    real(wp_), dimension(nlmt) :: chlm
+    integer :: nlm,ierr,npar
+
+    cst2=1.0_wp_-rbx
+    if(cst2<cst2min) cst2=0.0_wp_
+    call SpitzFuncCoeff(amu,Zeff,fc)
+    nlm=nlmt
+    call profil(0,tjp,njpt,tlm,nlmt,ch,ksp,ksp,rhop,nlm,chlm,ierr)
+    if(ierr>0) write(*,*) '    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
+  end subroutine setcdcoeff_ncl
+
+  subroutine eccdeff(yg,anpl,anprre,dens,amu,ex,ey,ez,nhmn,nhmx,       &
+        ithn,cst2,fcur,eccdpar,effjcd,iokhawa,ierr)
+    use const_and_precisions, only : pi,qesi=>e_,mesi=>me_,              &
+        vcsi=>c_,qe=>ecgs_,me=>mecgs_,vc=>ccgs_,mc2=>mc2_
+    use quadpack, only : dqagsmv
+    implicit none
+! 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_
+    real(wp_), parameter :: dumin=1.0e-6_wp_
+    integer, parameter :: lw=5000,liw=lw/4
+! arguments
+    integer :: i,nhmn,nhmx,ithn,iokhawa,ierr
+    real(wp_) :: yg,anpl,anprre,dens,amu,cst2,effjcd
+    real(wp_), dimension(:) :: eccdpar
+    complex(wp_) :: ex,ey,ez
+! local variables
+    integer :: nhn,neval,ier,last,npar
+    integer, dimension(liw) :: iw
+    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+size(eccdpar)) :: apar
+    real(wp_), dimension(0:1) :: uleft,uright
+! common/external functions/variables
+    real(wp_), external :: fcur
+!
+!     effjpl = <J_parallel>/<p_d> /(B_min/<B>)  [A m /W]
+!
+    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)
+
+    npar=size(apar)
+    apar(nfpp+1:npar) = eccdpar
+
+    anpl2=anpl*anpl
+
+    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
+
+      rdu2=anpl2+ygn2-1.0_wp_
+
+      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
+
+      uu1=uplm
+      uu2=uplp
+      xx1=amu*(anpl*uu1+ygn-1.0_wp_)
+      xx2=amu*(anpl*uu2+ygn-1.0_wp_)
+
+      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.le.dumin) cycle
+
+      apar(12) = dble(nhn)
+      apar(13) = ygn
+
+      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
+
+      rdu2t=cst2*anpl2+ygn2-1.0_wp_
+
+      if (rdu2t.gt.0.0_wp_.and.cst2.gt.0.0_wp_) then
+!
+!       resonance curve crosses the trapping region
+!
+        iokhawa=1
+        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
+!
+!       resonance curve does not cross the trapping region
+!
+        iokhawa=0
+        uleft(0)=uplm
+        uright(0)=uplp
+      end if
+
+      resj=0.0_wp_
+      do i=0,1
+        resji=0.0_wp_
+        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
+
+    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
+
+  end subroutine eccdeff
+
+  function fpp(upl,extrapar,npar)
+!
+!     computation of integral for power density, integrand function fpp
+!
+!     ith=0 : polarization term = const
+!     ith=1 : polarization term  Larmor radius expansion to lowest order
+!     ith=2 : full polarization term (J Bessel)
+!
+!     integration variable upl passed explicitly. other variables passed
+!     as array of extra parameters of length npar=size(extrapar)
+!     
+!      extrapar(1) = yg
+!      extrapar(2) = anpl
+!      extrapar(3) = amu
+!      extrapar(4) = Re(anprw)
+!      extrapar(5) = Re(ex)
+!      extrapar(6) = Im(ex)
+!      extrapar(7) = Re(ey)
+!      extrapar(8) = Im(ey)
+!      extrapar(9) = Re(ez)
+!      extrapar(10) = Im(ez)
+!      extrapar(11) = double(ithn)
+!      extrapar(12) = double(nhn)
+!      extrapar(13) = ygn
+!
+    use const_and_precisions, only : ui=>im
+    use math, only : fact
+    implicit none
+! arguments
+    integer :: npar
+    real(wp_) :: upl,fpp
+    real(wp_), dimension(npar) :: extrapar
+! local variables
+    integer :: ithn,nhn,nm,np
+    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
+
+    yg=extrapar(1)
+    anpl=extrapar(2)
+    amu=extrapar(3)
+    anprre=extrapar(4)
+    ex=cmplx(extrapar(5),extrapar(6),wp_)
+    ey=cmplx(extrapar(7),extrapar(8),wp_)
+    ez=cmplx(extrapar(9),extrapar(10),wp_)
+    ithn=int(extrapar(11))
+    nhn=int(extrapar(12))
+    ygn=extrapar(13)
+
+    gam=anpl*upl+ygn
+    upr2=gam*gam-1.0_wp_-upl*upl
+    ee=exp(-amu*(gam-1))
+
+!    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
+!     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)
+          thn2=(0.5_wp_*cth*abs(emxy+ez*anprre*upl/ygn))**2
+          thn2u=upr2*thn2
+        else
+!     Full polarization term
+          nm=nhn-1
+          np=nhn+1
+          ajbnm=dbesjn(nm, bb)
+          ajbnp=dbesjn(np, bb)
+          ajbn=dbesjn(nhn, bb)
+          thn2u=(abs(ez*ajbn*upl+upr*(ajbnp*epxy+ajbnm*emxy)/2.0_wp_))**2
+        end if
+      end if
+    end if
+
+    fpp=ee*thn2u
+  end function fpp
+
+  function fjch(upl,extrapar,npar)
+!
+!     computation of integral for current density
+!     integrand for Cohen model with trapping
+!
+!     integration variable upl passed explicitly. Other variables passed
+!     as array of extra parameters of length npar=size(extrapar).
+!     variables with index 1..nfpp must be passed to fpp
+!     variable with index nfpp+1 is zeff
+!     variables with index gt nfpp+1 are specific of the cd model
+!     
+!      extrapar(2) = anpl
+!      extrapar(4) = Re(anprw)
+!      extrapar(13) = ygn
+!
+!      extrapar(14) = zeff
+!      extrapar(15) = rb
+!      extrapar(16) = alams
+!      extrapar(17) = pa
+!      extrapar(18) = fp0s
+!
+    use conical, only : fconic
+    implicit none
+! arguments
+    integer :: npar
+    real(wp_) :: upl,fjch
+    real(wp_), dimension(npar) :: extrapar
+! local variables
+    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
+
+    anpl=extrapar(2)
+    anprre=extrapar(4)
+    ygn=extrapar(13)
+    zeff=extrapar(nfpp1)
+    rb=extrapar(nfpp2)
+    alams=extrapar(nfpp3)
+    pa=extrapar(nfpp4)
+    fp0s=extrapar(nfpp5)
+
+    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
+    xib=1.0_wp_-xi
+    xibi=1.0_wp_/xib
+    fu2b=1.0_wp_+xib*u2
+    fu2=1.0_wp_+xi*u2
+    gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2)
+    gg=u*gu/z5
+    dgg=(gu+u2*(2.0_wp_/fu2b**(1.0_wp_+xibi)/sqrt(fu2)-xi*gu/fu2))/z5
+
+    alam=sqrt(1.0_wp_-upr2/u2/rb)
+    fp0=fconic(alam,pa,0)
+    dfp0=-(pa*pa/2.0_wp_+0.125_wp_)
+    if (alam.lt.1.0_wp_) dfp0=-fconic(alam,pa,1)/sqrt(1.0_wp_-alam**2)
+    fh=alam*(1.0_wp_-alams*fp0/(alam*fp0s))
+    dfhl=1.0_wp_-alams*dfp0/fp0s
+
+    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(1:nfpp),nfpp)
+
+  end function fjch
+
+  function fjch0(upl,extrapar,npar)
+!
+!     computation of integral for current density
+!     integrand for Cohen model without trapping    
+!
+!     integration variable upl passed explicitly. Other variables passed
+!     as array of extra parameters of length npar=size(extrapar).
+!     variables with index 1..nfpp must be passed to fpp
+!     variable with index nfpp+1 is zeff
+!     variables with index gt nfpp+1 are specific of the cd model
+!     
+!      extrapar(2) = anpl
+!      extrapar(13) = ygn
+!
+!      extrapar(14) = zeff
+!
+    implicit none
+! arguments
+    real(wp_) :: upl,fjch0
+    integer :: npar
+    real(wp_), dimension(npar) :: extrapar
+! local variables
+    real(wp_) :: anpl,ygn,zeff,gam,u2,u,z5,xi,xib,xibi,fu2b,fu2,gu,gg,dgg,eta
+!
+    anpl=extrapar(2)
+    ygn=extrapar(13)
+    zeff=extrapar(nfpp1)
+
+    gam=anpl*upl+ygn
+    u2=gam*gam-1.0_wp_
+    u=sqrt(u2)
+    z5=Zeff+5.0_wp_
+    xi=1.0_wp_/z5**2
+    xib=1.0_wp_-xi
+    xibi=1.0_wp_/xib
+    fu2b=1.0_wp_+xib*u2
+    fu2=1.0_wp_+xi*u2
+    gu=(1.0_wp_-1.0_wp_/fu2b**xibi)/sqrt(fu2)
+    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(1:nfpp),nfpp)
+
+  end function fjch0
+
+  function fjncl(upl,extrapar,npar)
+!
+!     computation of integral for current density
+!     integrand for momentum conserv. model K(u) from Maruschenko
+!     gg=F(u)/u with F(u) as in Cohen paper
+!
+!     integration variable upl passed explicitly. Other variables passed
+!     as array of extra parameters of length npar=size(extrapar).
+!     variables with index 1..nfpp must be passed to fpp
+!     variable with index nfpp+1 is zeff
+!     variables with index gt nfpp+1 are specific of the cd model
+!     
+!      extrapar(2) = anpl
+!      extrapar(3) = amu
+!      extrapar(13) = ygn
+!
+!      extrapar(14) = zeff
+!      extrapar(15) = fc
+!      extrapar(16) = rbx
+!      extrapar(17:16+(npar-16)/2) = tlm
+!      extrapar(17+(npar-16)/2:npar) = chlm
+!
+    use dierckx, only : splev,splder
+    implicit none
+! arguments
+    integer :: npar
+    real(wp_) :: upl,fjncl
+    real(wp_), dimension(npar) :: extrapar
+! local variables
+    integer :: nlm
+    real(wp_) :: anpl,amu,ygn,zeff,fc,rbx,gam,u2,u,upr2,  &
+      bth,uth,fk,dfk,alam,fu,dfu,eta
+! local variables
+    integer :: ier
+    real(wp_), dimension((npar-nfpp3)/2) :: wrk
+    real(wp_), dimension(1) :: xs,ys
+!
+    anpl=extrapar(2)
+    amu=extrapar(3)
+    ygn=extrapar(13)
+    zeff=extrapar(nfpp1)
+    fc=extrapar(nfpp2)
+    rbx=extrapar(nfpp3)
+
+    gam=anpl*upl+ygn
+    u2=gam*gam-1.0_wp_
+    u=sqrt(u2)
+    upr2=u2-upl*upl
+    bth=sqrt(2.0_wp_/amu)
+    uth=u/bth
+    call GenSpitzFunc(Zeff,fc,uth,u,gam,fk,dfk)
+    fk=fk*(4.0_wp_/amu**2)
+    dfk=dfk*(2.0_wp_/amu)*bth
+      
+    alam=upr2/u2/rbx
+    xs(1)=alam
+    nlm=(npar-nfpp3)/2
+!
+!      extrapar(17:16+(npar-16)/2) = tlm
+!      extrapar(17+(npar-16)/2:npar) = chlm
+!
+    call splev(extrapar(nfpp4:nfpp3+nlm),nlm,extrapar(nfpp4+nlm:npar),3,     &
+               xs(1),ys(1),1,ier)
+    fu=ys(1)
+    call splder(extrapar(nfpp4:nfpp3+nlm),nlm,extrapar(nfpp4+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(1:nfpp),nfpp)
+  end function fjncl
+
+  SUBROUTINE GenSpitzFunc(Zeff,fc,u,q,gam, K,dKdu)
+!=======================================================================
+! Author:  N.B.Marushchenko
+! June 2005: as start point the subroutine of Ugo Gasparino (198?) 
+!            SpitzFunc() is taken and modified.
+!            1. adapted to the Fortran-95
+!            2. derivative of Spitzer function is added
+!            3. separation for 2 brunches is done: 
+!               1st is referenced as 'with conservation of the moment',
+!               2nd - as 'high speed limit'. 
+!               The last one is taken from the Lin-Liu formulation 
+!               (Phys.Plasmas 10 (2003) 4064) with K = F*fc.
+!               The asymptotical high speed limit (Taguchi-Fisch model)
+!               is also included as the reference case.
+! Feb. 2008: non-relativ. version is replaced by the relativistic one;
+!            the method is the the same, but the trial-function is 
+!            based on the relativistic formulation.
+!            The relativistic corrections for the collisional operator 
+!            up to the second order, i.e. (1/mu)**2, are applied.
+! Sep. 2008: generalized Spitzer function for arbitrary collisionality
+!            is implemented. The model is based on the concept of 
+!            the "effective trapped particles fraction". 
+!            The different.-integral kinetic equation for the generalized
+!            Spitzer function is produced with help of subroutines 
+!            ArbColl_TrappFract_Array and ArbColl_SpitzFunc_Array,
+!            where the subroutines of H. Maassberg are called).
+!========================================================================
+! Spitzer function with & w/o trapped particle effects is given by:
+!
+!    K(x) = x/gamma*(d1*x+d2*x^2+d4*x^3+d4*x^4),
+! 
+! where x = v/v_th and gamma=1 for non-relativistic version (Ugo),
+! or    x = p/p_th for relativistic version (N.M., February 2008).
+! Note, that somewhere the function F(x) instead of K(x) is applied,
+!
+!    F(x) = K(x)/fc.
+!
+! Numerical inversion of the 5x5 symmetric matrix obtained from the
+! generalized Spitzer problem (see paper of Taguchi for the equation
+! and paper of Hirshman for the variational approach bringing to the
+! matrix to be inverted).
+! 
+! The numerical method used is an improved elimination scheme
+! (Banachiewiczs-Cholesky-Crout method).
+! This method is particularly simple for symmetric matrix.
+! As a reference see "Mathematical Handbook" by Korn & Korn, p.635-636.
+! 
+! Refs.:  1. S.P. Hirshman, Phys. Fluids 23 (1980) 1238
+!         2. M. Rome' et al., Plasma Phys. Contr. Fus. 40 (1998) 511
+!         3. N.B. Marushchenko et al., Fusion Sci. Technol. 55 (2009) 180
+!========================================================================
+! INPUTS:
+!  u    - p/sqrt(2mT)
+!  q    - p/mc;
+!  gam  - relativistic factor;
+!  Zeff - effective charge;
+!  fc   - fraction of circulating particles.
+!
+! OUTPUTS:
+!   K   - Spitzer's function
+!  dKdu = dK/du, i.e. its derivative over normalized momentum
+!=======================================================================
+    use const_and_precisions, only : comp_eps
+    IMPLICIT NONE
+    REAL(wp_), INTENT(in)  :: Zeff,fc,u,q,gam
+    REAL(wp_), INTENT(out) :: K,dKdu
+    REAL(wp_) :: gam1,gam2,gam3
+
+     K   = 0
+    dKdu = 0
+    IF (u < comp_eps) RETURN
+
+    SELECT CASE(adj_appr(2))
+    CASE('m')  !--------------- momentum conservation ------------------!
+      gam1 = gam                                                        !
+      IF (adj_appr(4) == 'n') gam1 = 1                                  !
+      gam2 = gam1*gam1                                                  !
+      gam3 = gam1*gam2                                                  !
+       K   = u/gam1*u*(sfd(1)+u*(sfd(2)+u*(sfd(3)+u*sfd(4))))           !
+      dKdu = u/gam3*  (sfd(1)*(1+  gam2)+u*(sfd(2)*(1+2*gam2)+ &        !
+                    u*(sfd(3)*(1+3*gam2)+u* sfd(4)*(1+4*gam2))))        !
+    !--------------------- end momentum conservation -------------------!
+    CASE('h')  !---------------- high-speed-limit ----------------------!
+      IF (adj_appr(4) == 'n') THEN     !- non-relativ. asymptotic form -!
+         K   =   u**4 *fc/(Zeff+1+4*fc) !- (Taguchi-Fisch model)       -!
+        dKdu = 4*u**3 *fc/(Zeff+1+4*fc)                                 !
+      ELSEIF (adj_appr(4) == 'r') THEN !- relativistic, Lin-Liu form.  -!
+        CALL SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)          !
+      ENDIF                                                             !
+    CASE default   !----------------------------------------------------!
+      PRINT*,'GenSpitzFunc: WARNING! Spitzer function is not defined.'
+      RETURN
+    END SELECT
+  END SUBROUTINE GenSpitzFunc
+
+  SUBROUTINE SpitzFuncCoeff(mu,Zeff,fc)
+!=======================================================================
+! Calculates the matrix coefficients required for the subroutine 
+! "GenSpitzFunc", where the Spitzer function is defined through the
+! variational principle.
+! 
+! Weakly relativistic (upgraded) version (10.09.2008). 
+! Apart of the non-relativistic matrix coefficients, taken from the
+! old subroutine of Ugo Gasparino, the relativistic correction written
+! as series in 1/mu^n (mu=mc2/T) powers is added. Two orders are taken 
+! into account, i.e. n=0,1,2.
+!
+! In this version, the coefficients "oee", i.e. Omega_ij, are formulated
+! for arbitrary collisionality. 
+!
+! INPUT VARIABLES:
+! rho  = sqrt(SS) with SS - flux-surface label (norm. magn. flux)
+! ne   - density, 1/m^3
+! mu   - mc2/Te
+! Zeff - effective charge
+! fc   - fraction of circulating particles
+!
+! OUTPUT VARIABLES (defined as a global ones): 
+! sfd(1),...,sfd(4) - coefficients of the polynomial expansion of the 
+!                 "Spitzer"-function (the same as in the Hirshman paper)
+!=======================================================================
+    use const_and_precisions, only : mc2_
+    IMPLICIT NONE
+    REAL(wp_), INTENT(in) :: mu,Zeff,fc
+    INTEGER   :: n,i,j
+    REAL(wp_) :: rtc,rtc1,y,tn(1:nre)
+    REAL(wp_) :: m(0:4,0:4),g(0:4)
+    REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, &
+                       gam22,gam32,gam42,gam02, &
+                             gam33,gam43,gam03, &
+                                   gam44,gam04,gam00
+    REAL(wp_) :: alp12,alp13,alp14,alp10, &
+                       alp23,alp24,alp20, &
+                             alp34,alp30,alp40
+    REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0
+    LOGICAL   :: renew,rel,newmu,newZ,newfc
+    REAL(wp_), SAVE :: sfdx(1:4) = 0
+    REAL(wp_), SAVE :: mu_old =-1, Zeff_old =-1, fc_old =-1
+
+    rel   = mu < mc2_
+    newmu = abs(mu  -mu_old  ) > delta*mu
+    newZ  = abs(Zeff-Zeff_old) > delta*Zeff
+    newfc = abs(fc  -fc_old  ) > delta*fc
+    SELECT CASE(adj_appr(1))
+    CASE ('l','c')
+      renew = (newmu .and. rel) .OR. newZ .OR. newfc
+    END SELECT
+    IF (.not.renew) THEN
+      sfd(:) = sfdx(:)
+      RETURN
+    ENDIF
+
+    tn(:) = 0
+    IF (adj_appr(4) == 'r') THEN
+      IF (nre > 0) THEN
+        !mu = min(mu,1.e3*mc2_)
+        tn(1) = 1/mu
+        DO n=2,min(2,nre)
+          tn(n) = tn(n-1)/mu
+        ENDDO
+      ENDIF
+    ENDIF
+
+    SELECT CASE(adj_appr(1))
+    CASE ('l','c')     !---- both classical & collisionless limits  ----!
+      rtc = (1-fc)/fc;   rtc1 = rtc+1                                   !
+      !---                                                              !
+      DO i=0,4                                                          !
+        g(i) = vp_g(i,0)                                                !
+        DO n=1,min(2,nre)                                               !
+          g(i) = g(i) + tn(n)*vp_g(i,n)                                 !
+        ENDDO                                                           !
+        !---                                                            !
+        DO j=0,4                                                        !
+          IF (i == 0 .or. j == 0 .or. j >= i) THEN                      !
+              y =      vp_mee(i,j,0) + rtc *vp_oee(i,j,0) + &           !
+                                  Zeff*rtc1*vp_mei(i,j,0)               !
+            DO n=1,min(2,nre)                                           !
+              y = y + (vp_mee(i,j,n) + rtc *vp_oee(i,j,n) + &           !
+                                  Zeff*rtc1*vp_mei(i,j,n))*tn(n)        !
+            ENDDO                                                       !
+            m(i,j) = y                                                  !
+          ENDIF                                                         !
+        ENDDO                                                           !
+      ENDDO                                                             !
+      DO i=2,4                                                          !
+        DO j=1,i-1                                                      !
+          m(i,j) = m(j,i)                                               !
+        ENDDO                                                           !
+      ENDDO                                                             !
+      m(0,0) = 0                                                        !
+    CASE default       !------------------------------------------------!
+      PRINT*,'Green_Func: WARNING! Adjoint approach is not defined.'
+      RETURN
+    END SELECT
+
+    gam11 = m(1,1)
+    gam21 = m(2,1)
+    gam31 = m(3,1)
+    gam41 = m(4,1)
+    gam01 = m(0,1)
+
+    alp12 = m(1,2)/m(1,1)
+    alp13 = m(1,3)/m(1,1)
+    alp14 = m(1,4)/m(1,1)
+    alp10 = m(1,0)/m(1,1)
+
+    gam22 = m(2,2)-gam21*alp12
+    gam32 = m(3,2)-gam31*alp12
+    gam42 = m(4,2)-gam41*alp12
+    gam02 = m(0,2)-gam01*alp12
+
+    alp23 = gam32/gam22
+    alp24 = gam42/gam22
+    alp20 = gam02/gam22
+
+    gam33 = m(3,3)-gam31*alp13-gam32*alp23
+    gam43 = m(4,3)-gam41*alp13-gam42*alp23
+    gam03 = m(0,3)-gam01*alp13-gam02*alp23
+
+    alp34 = gam43/gam33
+    alp30 = gam03/gam33
+
+    gam44 = m(4,4)-gam41*alp14-gam42*alp24-gam43*alp34
+    gam04 = m(0,4)-gam01*alp14-gam02*alp24-gam03*alp34
+
+    alp40 = gam04/gam44
+
+    gam00 = m(0,0)-gam01*alp10-gam02*alp20-gam03*alp30-gam04*alp40
+
+    bet1 =  g(1)/m(1,1)
+    bet2 = (g(2)-gam21*bet1)/gam22
+    bet3 = (g(3)-gam31*bet1-gam32*bet2)/gam33
+    bet4 = (g(4)-gam41*bet1-gam42*bet2-gam43*bet3)/gam44
+    bet0 = (g(0)-gam01*bet1-gam02*bet2-gam03*bet3-gam04*bet4)/gam00
+
+    d0     = bet0
+    sfd(4) = bet4-alp40*d0
+    sfd(3) = bet3-alp30*d0-alp34*sfd(4)
+    sfd(2) = bet2-alp20*d0-alp24*sfd(4)-alp23*sfd(3)
+    sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2)
+
+    fc_old   = fc
+    mu_old   = mu
+    Zeff_old = Zeff
+    sfdx(1:4) = sfd(1:4)
+
+  END SUBROUTINE SpitzFuncCoeff
+
+  SUBROUTINE SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)
+!=======================================================================
+! Calculates the "Spitzer function" in high velocity limit, relativistic 
+! formulation: Lin-Liu et al., Phys.Pl. (2003),v10, 4064, Eq.(33).
+! 
+! Inputs:
+! Zeff - effective charge
+! fc   - fraction of circulating electrons
+! u    - p/(m*vte)
+! q    - p/mc
+! gam  - relativ. factor
+! 
+! Outputs:
+!  K   - Spitzer function
+! dKdu - its derivative
+!=======================================================================
+    use const_and_precisions, only : zero,one
+    use numint, only : quanc8
+    IMPLICIT NONE
+    REAL(wp_), INTENT(in)  :: Zeff,fc,u,q,gam
+    REAL(wp_), INTENT(out) :: K,dKdu
+    INTEGER :: nfun
+    REAL(wp_) :: gam2,err,flag,Integr
+    REAL(wp_), PARAMETER :: a = zero, b = one, rtol = 1e-4_wp_, atol = 1e-12_wp_
+
+    r2 = (1+Zeff)/fc   ! global parameter needed for integrand, HSL_f(t)
+
+    IF (u < 1e-2) THEN
+       K   =   u**4/(r2+4)
+      dKdu = 4*u**3/(r2+4)
+      RETURN
+    ENDIF
+
+    q2  = q*q       ! for the integrand, HSL_f
+    gp1 = gam+1     ! ..
+    CALL quanc8(HSL_f,zero,one,atol,rtol,Integr,err,nfun,flag)
+
+    gam2 = gam*gam
+    K   = u**4 * Integr
+    dKdu = (u/gam)**3 * (1-r2*gam2*Integr)
+  END SUBROUTINE SpitzFunc_HighSpeedLimit
+
+  FUNCTION HSL_f(t) RESULT(f)
+!=======================================================================
+! Integrand for the high-speed limit approach (Lin-Liu's formulation)
+!=======================================================================
+    IMPLICIT NONE
+    REAL(wp_), INTENT(in) :: t
+    REAL(wp_) :: f,g
+    g = sqrt(1+t*t*q2)
+    f = t**(3+r2)/g**3 * (gp1/(g+1))**r2
+  END FUNCTION HSL_f
+
+end module eccd

src/eierf.f90

@@ -0,0 +1,906 @@
+module eierf
+
+  use const_and_precisions, only : wp_, zero, one
+  implicit none
+  real(wp_), parameter, private :: half=0.5_wp_, two=2.0_wp_, three=3.0_wp_, &
+        four=4.0_wp_, six=6.0_wp_, twelve=12._wp_, sixten=16.0_wp_, &
+        two4=24.0_wp_, fourty=40.0_wp_
+
+contains
+
+! ======================================================================
+! nist guide to available math software.
+! fullsource for module ei from package specfun.
+! retrieved from netlib on fri mar 26 05:52:39 1999.
+! ======================================================================
+  subroutine calcei(arg,result,intt)
+!----------------------------------------------------------------------
+!
+! this fortran 77 packet computes the exponential integrals ei(x),
+!  e1(x), and  exp(-x)*ei(x)  for real arguments  x  where
+!
+!       integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
+!  ei(x) =
+!      -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
+!
+!  and where the first integral is a principal value integral.
+!  the packet contains three function type subprograms: ei, eone,
+!  and expei;  and one subroutine type subprogram: calcei.  the
+!  calling statements for the primary entries are
+!
+!             y = ei(x),            where  x /= 0,
+!
+!             y = eone(x),          where  x > 0,
+!  and
+!             y = expei(x),         where  x /= 0,
+!
+!  and where the entry points correspond to the functions ei(x),
+!  e1(x), and exp(-x)*ei(x), respectively.  the routine calcei
+!  is intended for internal packet use only, all computations within
+!  the packet being concentrated in this routine.  the function
+!  subprograms invoke calcei with the fortran statement
+!     call calcei(arg,result,intt)
+!  where the parameter usage is as follows
+!
+! function                  parameters for calcei
+!   call                 arg             result         intt
+!
+!  ei(x)              x /= 0          ei(x)            1
+!  eone(x)            x > 0          -ei(-x)           2
+!  expei(x)           x /= 0          exp(-x)*ei(x)    3
+!----------------------------------------------------------------------
+    implicit none
+    integer, intent(in) :: intt
+    real(wp_), intent(in) :: arg
+    real(wp_), intent(out) :: result
+    integer :: i
+    real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq
+    real(wp_), dimension(10) :: px,qx
+!----------------------------------------------------------------------
+!  mathematical constants
+!   exp40 = exp(40)
+!   x0 = zero of ei
+!   x01/x11 + x02 = zero of ei to extra precision
+!----------------------------------------------------------------------
+    real(wp_), parameter :: p037=0.037_wp_, &
+         exp40=2.3538526683701998541e17_wp_, x01=381.5_wp_, x11=1024.0_wp_, &
+         x02=-5.1182968633365538008e-5_wp_, x0=3.7250741078136663466e-1_wp_
+!----------------------------------------------------------------------
+! machine-dependent constants
+!----------------------------------------------------------------------
+    real(wp_), parameter :: xinf=1.79e+308_wp_,xmax=716.351_wp_,xbig=701.84_wp_
+!----------------------------------------------------------------------
+! coefficients  for -1.0 <= x < 0.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(7), parameter :: &
+       a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, &
+           1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, &
+           1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, &
+           5.0196785185439843791020_wp_/)
+    real(wp_), dimension(6), parameter :: &
+       b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, &
+           8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, &
+           1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/)
+!----------------------------------------------------------------------
+! coefficients for -4.0 <= x < -1.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(9), parameter :: &
+       c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, &
+           7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, &
+           1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, &
+           1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, &
+           1.737331760720576030932e-8_wp_/), &
+       d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, &
+           4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, &
+           2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, &
+           9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, &
+           1.000000000000000000000e+0_wp_/)
+!----------------------------------------------------------------------
+! coefficients for x < -4.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, &
+           1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, &
+           1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, &
+           1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, &
+           5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/),&
+       f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, &
+           5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, &
+           2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, &
+           1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, &
+           5.4199632588522559414924e+1_wp_,1.0_wp_/)
+!----------------------------------------------------------------------
+!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
+!----------------------------------------------------------------------
+    real(wp_), dimension(4), parameter :: &
+       plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, &
+             -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), &
+       qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, &
+             -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/)
+!----------------------------------------------------------------------
+! coefficients for  0.0 < x < 6.0,
+!  ratio of chebyshev polynomials
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, &
+           -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, &
+           -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, &
+            3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, &
+            1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),&
+       q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, &
+            1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, &
+            6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, &
+            5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, &
+            9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for 6.0 <= x < 12.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, &
+           -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, &
+            1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, &
+            1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, &
+            9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, &
+            3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, &
+            7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, &
+            4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, &
+            1.249884822712447891440_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for 12.0 <= x < 24.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, &
+            -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, &
+            -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, &
+             2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, &
+            -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, &
+             5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, &
+             4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, &
+             1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, &
+             1.001533852045342697818_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for  x >= 24.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, &
+            -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, &
+            -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, &
+            -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, &
+            -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, &
+             1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, &
+             7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, &
+            -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, &
+             1.99999999999048104167_wp_/)
+!----------------------------------------------------------------------
+    x = arg
+    if (x == zero) then
+      ei = -xinf
+      if (intt == 2) ei = -ei
+    else if ((x < zero) .or. (intt == 2)) then
+!----------------------------------------------------------------------
+! calculate ei for negative argument or for e1.
+!----------------------------------------------------------------------
+      y = abs(x)
+      if (y <= one) then
+        sump = a(7) * y + a(1)
+        sumq = y + b(1)
+        do i = 2, 6
+          sump = sump * y + a(i)
+          sumq = sumq * y + b(i)
+        end do
+        ei = log(y) - sump / sumq
+        if (intt == 3) ei = ei * exp(y)
+      else if (y <= four) then
+        w = one / y
+        sump = c(1)
+        sumq = d(1)
+        do i = 2, 9
+          sump = sump * w + c(i)
+          sumq = sumq * w + d(i)
+        end do
+        ei = - sump / sumq
+        if (intt /= 3) ei = ei * exp(-y)
+      else
+        if ((y > xbig) .and. (intt < 3)) then
+          ei = zero
+        else
+          w = one / y
+          sump = e(1)
+          sumq = f(1)
+          do i = 2, 10
+            sump = sump * w + e(i)
+            sumq = sumq * w + f(i)
+          end do
+          ei = -w * (one - w * sump / sumq )
+          if (intt /= 3) ei = ei * exp(-y)
+        end if
+      end if
+      if (intt == 2) ei = -ei
+    else if (x < six) then
+!----------------------------------------------------------------------
+!  to improve conditioning, rational approximations are expressed
+!    in terms of chebyshev polynomials for 0 <= x < 6, and in
+!    continued fraction form for larger x.
+!----------------------------------------------------------------------
+      t = x + x
+      t = t / three - two
+      px(1) = zero
+      qx(1) = zero
+      px(2) = p(1)
+      qx(2) = q(1)
+      do i = 2, 9
+        px(i+1) = t * px(i) - px(i-1) + p(i)
+        qx(i+1) = t * qx(i) - qx(i-1) + q(i)
+      end do
+      sump = half * t * px(10) - px(9) + p(10)
+      sumq = half * t * qx(10) - qx(9) + q(10)
+      frac = sump / sumq
+      xmx0 = (x - x01/x11) - x02
+      if (abs(xmx0) >= p037) then
+        ei = log(x/x0) + xmx0 * frac
+        if (intt == 3) ei = exp(-x) * ei
+      else
+!----------------------------------------------------------------------
+! special approximation to  ln(x/x0)  for x close to x0
+!----------------------------------------------------------------------
+        y = xmx0 / (x + x0)
+        ysq = y*y
+        sump = plg(1)
+        sumq = ysq + qlg(1)
+        do i = 2, 4
+          sump = sump*ysq + plg(i)
+          sumq = sumq*ysq + qlg(i)
+        end do
+        ei = (sump / (sumq*(x+x0)) + frac) * xmx0
+        if (intt == 3) ei = exp(-x) * ei
+      end if
+    else if (x < twelve) then
+      frac = zero
+      do i = 1, 9
+        frac = s(i) / (r(i) + x + frac)
+      end do
+      ei = (r(10) + frac) / x
+      if (intt /= 3) ei = ei * exp(x)
+    else if (x <= two4) then
+      frac = zero
+      do i = 1, 9
+        frac = q1(i) / (p1(i) + x + frac)
+      end do
+      ei = (p1(10) + frac) / x
+      if (intt /= 3) ei = ei * exp(x)
+    else
+      if ((x >= xmax) .and. (intt < 3)) then
+        ei = xinf
+      else
+        y = one / x
+        frac = zero
+        do i = 1, 9
+          frac = q2(i) / (p2(i) + x + frac)
+        end do
+        frac = p2(10) + frac
+        ei = y + y * y * frac
+        if (intt /= 3) then
+          if (x <= xmax-two4) then
+            ei = ei * exp(x)
+          else
+!----------------------------------------------------------------------
+! calculation reformulated to avoid premature overflow
+!----------------------------------------------------------------------
+            ei = (ei * exp(x-fourty)) * exp40
+          end if
+        end if
+      end if
+    end if
+    result = ei
+  end subroutine calcei
+
+  function ei(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   exponential integral  ei(x), where  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+    implicit none
+    integer :: intt
+    real(wp_) :: ei
+    real(wp_), intent(in) :: x
+    real(wp_) :: result
+!--------------------------------------------------------------------
+    intt = 1
+    call calcei(x,result,intt)
+    ei = result
+  end function ei
+
+  function expei(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   function  exp(-x) * ei(x), where  ei(x)  is the exponential
+!   integral, and  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+    implicit none
+    integer :: intt
+    real(wp_) :: expei
+    real(wp_), intent(in) :: x
+    real(wp_) :: result
+!--------------------------------------------------------------------
+    intt = 3
+    call calcei(x,result,intt)
+    expei = result
+  end function expei
+
+  function eone(x)
+!--------------------------------------------------------------------
+!
+! this function program computes approximate values for the
+!   exponential integral e1(x), where  x  is real.
+!
+!  author: w. j. cody
+!
+!  latest modification: january 12, 1988
+!
+!--------------------------------------------------------------------
+    implicit none
+    integer :: intt
+    real(wp_) :: eone
+    real(wp_), intent(in) :: x
+    real(wp_) :: result
+!--------------------------------------------------------------------
+    intt = 2
+    call calcei(x,result,intt)
+    eone = result
+  end function eone
+
+! ======================================================================
+! calcei3 = calcei for int=3
+! ======================================================================
+  subroutine calcei3(arg,result)
+!----------------------------------------------------------------------
+!
+! this fortran 77 packet computes the exponential integrals ei(x),
+!  e1(x), and  exp(-x)*ei(x)  for real arguments  x  where
+!
+!       integral (from t=-infinity to t=x) (exp(t)/t),  x > 0,
+!  ei(x) =
+!      -integral (from t=-x to t=infinity) (exp(t)/t),  x < 0,
+!
+!  and where the first integral is a principal value integral.
+!  the packet contains three function type subprograms: ei, eone,
+!  and expei;  and one subroutine type subprogram: calcei.  the
+!  calling statements for the primary entries are
+!
+!             y = ei(x),            where  x /= 0,
+!
+!             y = eone(x),          where  x > 0,
+!  and
+!             y = expei(x),         where  x /= 0,
+!
+!  and where the entry points correspond to the functions ei(x),
+!  e1(x), and exp(-x)*ei(x), respectively.  the routine calcei
+!  is intended for internal packet use only, all computations within
+!  the packet being concentrated in this routine.  the function
+!  subprograms invoke calcei with the fortran statement
+!     call calcei(arg,result,int)
+!  where the parameter usage is as follows
+!
+! function                  parameters for calcei
+!   call                 arg             result         int
+!
+!  ei(x)              x /= 0          ei(x)            1
+!  eone(x)            x > 0         -ei(-x)           2
+!  expei(x)           x /= 0          exp(-x)*ei(x)    3
+!----------------------------------------------------------------------
+    implicit none
+    real(wp_), intent(in) :: arg
+    real(wp_), intent(out) :: result
+    integer :: i
+    real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq
+    real(wp_), dimension(10) :: px,qx
+!----------------------------------------------------------------------
+!  mathematical constants
+!   exp40 = exp(40)
+!   x0 = zero of ei
+!   x01/x11 + x02 = zero of ei to extra precision
+!----------------------------------------------------------------------
+    real(wp_), parameter :: p037=0.037_wp_, &
+         x01=381.5_wp_, x11=1024.0_wp_, x02=-5.1182968633365538008e-5_wp_, &
+         x0=3.7250741078136663466e-1_wp_
+!----------------------------------------------------------------------
+! machine-dependent constants
+!----------------------------------------------------------------------
+    real(wp_), parameter :: xinf=1.79e+308_wp_
+!----------------------------------------------------------------------
+! coefficients  for -1.0 <= x < 0.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(7), parameter :: &
+       a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, &
+           1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, &
+           1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, &
+           5.0196785185439843791020_wp_/)
+    real(wp_), dimension(6), parameter :: &
+       b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, &
+           8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, &
+           1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/)
+!----------------------------------------------------------------------
+! coefficients for -4.0 <= x < -1.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(9), parameter :: &
+       c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, &
+           7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, &
+           1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, &
+           1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, &
+           1.737331760720576030932e-8_wp_/), &
+       d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, &
+           4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, &
+           2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, &
+           9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, &
+           1.000000000000000000000e+0_wp_/)
+!----------------------------------------------------------------------
+! coefficients for x < -4.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, &
+           1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, &
+           1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, &
+           1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, &
+           5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/), &
+       f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, &
+           5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, &
+           2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, &
+           1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, &
+           5.4199632588522559414924e+1_wp_,1.0_wp_/)
+!----------------------------------------------------------------------
+!  coefficients for rational approximation to ln(x/a), |1-x/a| < .1
+!----------------------------------------------------------------------
+    real(wp_), dimension(4), parameter :: &
+       plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, &
+             -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), &
+       qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, &
+             -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/)
+!----------------------------------------------------------------------
+! coefficients for  0.0 < x < 6.0,
+!  ratio of chebyshev polynomials
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, &
+           -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, &
+           -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, &
+            3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, &
+            1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),&
+       q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, &
+            1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, &
+            6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, &
+            5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, &
+            9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for 6.0 <= x < 12.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, &
+           -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, &
+            1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, &
+            1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, &
+            9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, &
+            3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, &
+            7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, &
+            4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, &
+            1.249884822712447891440_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for 12.0 <= x < 24.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, &
+            -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, &
+            -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, &
+             2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, &
+            -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, &
+             5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, &
+             4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, &
+             1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, &
+             1.001533852045342697818_wp_/)
+!----------------------------------------------------------------------
+! j-fraction coefficients for  x >= 24.0
+!----------------------------------------------------------------------
+    real(wp_), dimension(10), parameter :: &
+       p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, &
+            -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, &
+            -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, &
+            -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, &
+            -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/)
+    real(wp_), dimension(9), parameter :: &
+       q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, &
+             1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, &
+             7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, &
+            -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, &
+             1.99999999999048104167_wp_/)
+!----------------------------------------------------------------------
+    x = arg
+    if (x == zero) then
+      ei = -xinf
+    else if ((x < zero)) then
+!----------------------------------------------------------------------
+! calculate ei for negative argument or for e1.
+!----------------------------------------------------------------------
+      y = abs(x)
+      if (y <= one) then
+        sump = a(7) * y + a(1)
+        sumq = y + b(1)
+        do i = 2, 6
+          sump = sump * y + a(i)
+          sumq = sumq * y + b(i)
+        end do
+        ei = (log(y) - sump / sumq ) * exp(y)
+      else if (y <= four) then
+        w = one / y
+        sump = c(1)
+        sumq = d(1)
+        do i = 2, 9
+          sump = sump * w + c(i)
+          sumq = sumq * w + d(i)
+        end do
+        ei = - sump / sumq
+      else
+        w = one / y
+        sump = e(1)
+        sumq = f(1)
+        do i = 2, 10
+          sump = sump * w + e(i)
+          sumq = sumq * w + f(i)
+        end do
+        ei = -w * (one - w * sump / sumq )
+      end if
+    else if (x < six) then
+!----------------------------------------------------------------------
+!  to improve conditioning, rational approximations are expressed
+!    in terms of chebyshev polynomials for 0 <= x < 6, and in
+!    continued fraction form for larger x.
+!----------------------------------------------------------------------
+      t = x + x
+      t = t / three - two
+      px(1) = zero
+      qx(1) = zero
+      px(2) = p(1)
+      qx(2) = q(1)
+      do i = 2, 9
+        px(i+1) = t * px(i) - px(i-1) + p(i)
+        qx(i+1) = t * qx(i) - qx(i-1) + q(i)
+      end do
+      sump = half * t * px(10) - px(9) + p(10)
+      sumq = half * t * qx(10) - qx(9) + q(10)
+      frac = sump / sumq
+      xmx0 = (x - x01/x11) - x02
+      if (abs(xmx0) >= p037) then
+        ei = exp(-x) * ( log(x/x0) + xmx0 * frac )
+      else
+!----------------------------------------------------------------------
+! special approximation to  ln(x/x0)  for x close to x0
+!----------------------------------------------------------------------
+        y = xmx0 / (x + x0)
+        ysq = y*y
+        sump = plg(1)
+        sumq = ysq + qlg(1)
+        do i = 2, 4
+          sump = sump*ysq + plg(i)
+          sumq = sumq*ysq + qlg(i)
+        end do
+        ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0
+      end if
+    else if (x < twelve) then
+      frac = zero
+      do i = 1, 9
+        frac = s(i) / (r(i) + x + frac)
+      end do
+      ei = (r(10) + frac) / x
+    else if (x <= two4) then
+      frac = zero
+      do i = 1, 9
+        frac = q1(i) / (p1(i) + x + frac)
+      end do
+      ei = (p1(10) + frac) / x
+    else
+      y = one / x
+      frac = zero
+      do i = 1, 9
+        frac = q2(i) / (p2(i) + x + frac)
+      end do
+      frac = p2(10) + frac
+      ei = y + y * y * frac
+    end if
+    result = ei
+  end subroutine calcei3
+
+!  subroutine calerf(arg,result,jintt)
+!!------------------------------------------------------------------
+!!
+!! this packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x)
+!!   for a real argument  x.  it contains three function type
+!!   subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx),
+!!   and one subroutine type subprogram, calerf.  the calling
+!!   statements for the primary entries are:
+!!
+!!               y=erf(x)     (or   y=derf(x)),
+!!
+!!               y=erfc(x)    (or   y=derfc(x)),
+!!   and
+!!               y=erfcx(x)   (or   y=derfcx(x)).
+!!
+!!   the routine  calerf  is intended for internal packet use only,
+!!   all computations within the packet being concentrated in this
+!!   routine.  the function subprograms invoke  calerf  with the
+!!   statement
+!!
+!!      call calerf(arg,result,jintt)
+!!
+!!   where the parameter usage is as follows
+!!
+!!  function                     parameters for calerf
+!!   call              arg                  result          jintt
+!!
+!! erf(arg)      any real argument         erf(arg)          0
+!! erfc(arg)     abs(arg) < xbig           erfc(arg)         1
+!! erfcx(arg)    xneg < arg < xmax         erfcx(arg)        2
+!!
+!!*******************************************************************
+!!*******************************************************************
+!!
+!! Explanation of machine-dependent constants
+!!
+!!   XMIN   = the smallest positive floating-point number.
+!!   XINF   = the largest positive finite floating-point number.
+!!   XNEG   = the largest negative argument acceptable to ERFCX;
+!!            the negative of the solution to the equation
+!!            2*exp(x*x) = XINF.
+!!   XSMALL = argument below which erf(x) may be represented by
+!!            2*x/sqrt(pi)  and above which  x*x  will not underflow.
+!!            A conservative value is the largest machine number X
+!!            such that   1.0 + X = 1.0   to machine precision.
+!!   XBIG   = largest argument acceptable to ERFC;  solution to
+!!            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where
+!!            W(x) = exp(-x*x)/[x*sqrt(pi)].
+!!   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to
+!!            machine precision.  A conservative value is
+!!            1/[2*sqrt(XSMALL)]
+!!   XMAX   = largest acceptable argument to ERFCX; the minimum
+!!            of XINF and 1/[sqrt(pi)*XMIN].
+!!
+!!*******************************************************************
+!!*******************************************************************
+!!
+!! error returns
+!!
+!!  the program returns  erfc = 0      for  arg >= xbig;
+!!
+!!                   erfcx = xinf  for  arg < xneg;
+!!  and
+!!                   erfcx = 0     for  arg >= xmax.
+!!
+!!
+!! intrinsic functions required are:
+!!
+!! abs, aint, exp
+!!
+!!
+!!  author: w. j. cody
+!!      mathematics and computer science division
+!!      argonne national laboratory
+!!      argonne, il 60439
+!!
+!!  latest modification: march 19, 1990
+!!
+!!------------------------------------------------------------------
+!    implicit none
+!    real(wp_), intent(in) :: arg
+!    real(wp_), intent(out) :: result
+!    integer, intent(in) :: jintt
+!    integer :: i
+!    real(wp_) :: del,x,xden,xnum,y,ysq
+!!------------------------------------------------------------------
+!!  mathematical constants
+!!------------------------------------------------------------------
+!    real(wp_), parameter :: sqrpi=5.6418958354775628695e-1_wp_, &
+!       thresh=0.46875_wp_
+!!------------------------------------------------------------------
+!!  machine-dependent constants
+!!------------------------------------------------------------------
+!    real(wp_), parameter :: xinf=1.79e308_wp_,    &  ! ~huge
+!                            xneg=-26.628_wp_,     &  ! ?
+!                            xsmall=1.11e-16_wp_,  &  ! ~epsilon/2
+!                            xbig=26.543_wp_,      &  ! ?
+!                            xhuge=6.71e7_wp_,     &  ! ~1/sqrt(epsilon)
+!                            xmax=2.53e307_wp_        ! ?
+!!------------------------------------------------------------------
+!!  coefficients for approximation to  erf  in first interval
+!!------------------------------------------------------------------
+!    real(wp_), dimension(5), parameter :: &
+!       a=(/3.16112374387056560_wp_,1.13864154151050156e02_wp_, &
+!           3.77485237685302021e02_wp_,3.20937758913846947e03_wp_, &
+!           1.85777706184603153e-1_wp_/)
+!    real(wp_), dimension(4), parameter :: &
+!       b=(/2.36012909523441209e01_wp_,2.44024637934444173e02_wp_, &
+!           1.28261652607737228e03_wp_,2.84423683343917062e03_wp_/)
+!!------------------------------------------------------------------
+!!  coefficients for approximation to  erfc  in second interval
+!!------------------------------------------------------------------
+!    real(wp_), dimension(9), parameter :: &
+!       c=(/5.64188496988670089e-1_wp_,8.88314979438837594_wp_, &
+!           6.61191906371416295e01_wp_,2.98635138197400131e02_wp_, &
+!           8.81952221241769090e02_wp_,1.71204761263407058e03_wp_, &
+!           2.05107837782607147e03_wp_,1.23033935479799725e03_wp_, &
+!           2.15311535474403846e-8_wp_/)
+!    real(wp_), dimension(8), parameter :: &
+!       d=(/1.57449261107098347e01_wp_,1.17693950891312499e02_wp_, &
+!           5.37181101862009858e02_wp_,1.62138957456669019e03_wp_, &
+!           3.29079923573345963e03_wp_,4.36261909014324716e03_wp_, &
+!           3.43936767414372164e03_wp_,1.23033935480374942e03_wp_/)
+!!------------------------------------------------------------------
+!!  coefficients for approximation to  erfc  in third interval
+!!------------------------------------------------------------------
+!    real(wp_), dimension(6), parameter :: &
+!       p=(/3.05326634961232344e-1_wp_,3.60344899949804439e-1_wp_, &
+!           1.25781726111229246e-1_wp_,1.60837851487422766e-2_wp_, &
+!           6.58749161529837803e-4_wp_,1.63153871373020978e-2_wp_/)
+!    real(wp_), dimension(5), parameter :: &
+!       q=(/2.56852019228982242_wp_,1.87295284992346047_wp_, &
+!           5.27905102951428412e-1_wp_,6.05183413124413191e-2_wp_, &
+!           2.33520497626869185e-3_wp_/)
+!!------------------------------------------------------------------
+!    x = arg
+!    y = abs(x)
+!    if (y <= thresh) then
+!!------------------------------------------------------------------
+!!  evaluate  erf  for  |x| <= 0.46875
+!!------------------------------------------------------------------
+!      ysq = zero
+!      if (y > xsmall) ysq = y * y
+!      xnum = a(5)*ysq
+!      xden = ysq
+!      do i = 1, 3
+!        xnum = (xnum + a(i)) * ysq
+!        xden = (xden + b(i)) * ysq
+!      end do
+!      result = x * (xnum + a(4)) / (xden + b(4))
+!      if (jintt /= 0) result = one - result
+!      if (jintt == 2) result = exp(ysq) * result
+!      return
+!!------------------------------------------------------------------
+!!  evaluate  erfc  for 0.46875 <= |x| <= 4.0
+!!------------------------------------------------------------------
+!    else if (y <= four) then
+!      xnum = c(9)*y
+!      xden = y
+!      do i = 1, 7
+!        xnum = (xnum + c(i)) * y
+!        xden = (xden + d(i)) * y
+!      end do
+!      result = (xnum + c(8)) / (xden + d(8))
+!      if (jintt /= 2) then
+!        ysq = aint(y*sixten)/sixten
+!        del = (y-ysq)*(y+ysq)
+!        result = exp(-ysq*ysq) * exp(-del) * result
+!      end if
+!!------------------------------------------------------------------
+!!  evaluate  erfc  for |x| > 4.0
+!!------------------------------------------------------------------
+!    else if (y < xbig .or. (y < xmax .and. jintt == 2)) then
+!      ysq = one / (y * y)
+!      xnum = p(6)*ysq
+!      xden = ysq
+!      do i = 1, 4
+!        xnum = (xnum + p(i)) * ysq
+!        xden = (xden + q(i)) * ysq
+!      end do
+!      result = ysq *(xnum + p(5)) / (xden + q(5))
+!      result = (sqrpi -  result) / y
+!      if (jintt /= 2) then
+!        ysq = aint(y*sixten)/sixten
+!        del = (y-ysq)*(y+ysq)
+!        result = exp(-ysq*ysq) * exp(-del) * result
+!      end if
+!    else if (y >= xhuge) then
+!      result = sqrpi / y
+!    else
+!      result = zero
+!    end if
+!!------------------------------------------------------------------
+!!  fix up for negative argument, erf, etc.
+!!------------------------------------------------------------------
+!    if (jintt == 0) then
+!      result = (half - result) + half
+!      if (x < zero) result = -result
+!    else if (jintt == 1) then
+!      if (x < zero) result = two - result
+!    else
+!      if (x < zero) then
+!        if (x < xneg) then
+!          result = xinf
+!        else
+!          ysq = aint(x*sixten)/sixten
+!          del = (x-ysq)*(x+ysq)
+!          y = exp(ysq*ysq) * exp(del)
+!          result = (y+y) - result
+!        end if
+!      end if
+!    end if
+!  end subroutine calerf
+!
+!  function derf(x)
+!!--------------------------------------------------------------------
+!!
+!! this subprogram computes approximate values for erf(x).
+!!   (see comments heading calerf).
+!!
+!!   author/date: w. j. cody, january 8, 1985
+!!
+!!--------------------------------------------------------------------
+!    implicit none
+!    real(wp_) :: derf
+!    real(wp_), intent(in) :: x
+!    integer :: jintt
+!    real(wp_) :: result
+!!------------------------------------------------------------------
+!    jintt = 0
+!    call calerf(x,result,jintt)
+!    derf = result
+!  end function derf
+!
+!  function derfc(x)
+!!--------------------------------------------------------------------
+!!
+!! this subprogram computes approximate values for erfc(x).
+!!   (see comments heading calerf).
+!!
+!!   author/date: w. j. cody, january 8, 1985
+!!
+!!--------------------------------------------------------------------
+!    implicit none
+!    real(wp_) :: derfc
+!    real(wp_), intent(in) :: x
+!    integer :: jintt
+!    real(wp_) :: result
+!!------------------------------------------------------------------
+!    jintt = 1
+!    call calerf(x,result,jintt)
+!    derfc = result
+!  end function derfc
+!
+!  function derfcx(x)
+!!------------------------------------------------------------------
+!!
+!! this subprogram computes approximate values for exp(x*x) * erfc(x).
+!!   (see comments heading calerf).
+!!
+!!   author/date: w. j. cody, march 30, 1987
+!!
+!!------------------------------------------------------------------
+!    implicit none
+!    real(wp_) :: derfcx
+!    real(wp_), intent(in) :: x
+!    integer :: jintt
+!    real(wp_) :: result
+!!------------------------------------------------------------------
+!    jintt = 2
+!    call calerf(x,result,jintt)
+!    derfcx = result
+!  end function derfcx
+
+end module eierf

src/equilibrium.f90

@@ -0,0 +1,979 @@
+module equilibrium
+  use const_and_precisions, only : wp_
+  implicit none
+  
+  REAL(wp_), SAVE :: btaxis,rmaxis,zmaxis
+  REAL(wp_), SAVE :: btrcen ! used only for Jcd_ASTRA def.
+  REAL(wp_), SAVE :: rcen   ! computed as fpol(a)/btrcen
+  REAL(wp_), SAVE :: rmnm,rmxm,zmnm,zmxm
+  REAL(wp_), SAVE :: zbinf,zbsup
+  REAL(wp_), SAVE :: rup,zup,rlw,zlw
+
+  INTEGER, PARAMETER :: kspl=3,ksplp=kspl+1
+  ! === 2D spline psi(R,z), normalization and derivatives ==========
+  INTEGER, SAVE :: nsr, nsz
+  REAL(wp_), SAVE :: psia, psiant, psinop
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tr,tz
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: cceq,   cceq01, cceq10, &
+                                                cceq20, cceq02, cceq11
+  ! === 1D spline Fpol(psi) ========================================
+  ! INTEGER, SAVE :: npsiest
+  INTEGER, SAVE :: nsf
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tfp, cfp
+  REAL(wp_), SAVE :: fpolas
+  ! === 1D spline rhot(rhop), rhop(rhot), q(psi) ===================
+  ! computed on psinr,rhopnr [,rhotnr] arrays
+  INTEGER, SAVE :: nq,nrho
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psinr,rhopr,rhotr
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cq,crhop,crhot
+  REAL(wp_), SAVE :: phitedge,aminor
+  REAL(wp_), SAVE :: q0,qa,alq
+
+contains
+
+  subroutine read_eqdsk(filenm,rv,zv,psin,psia,psinr,fpol,q,rvac,rax,zax, &
+                        rbnd,zbnd,rlim,zlim,ipsinorm,idesc,ifreefmt,unit)
+    use const_and_precisions, only : one
+    use utils, only : get_free_unit
+    implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    real(wp_), intent(out) :: psia,rvac,rax,zax
+    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,psinr,fpol,q
+    real(wp_), dimension(:), allocatable, intent(out) :: rbnd,zbnd,rlim,zlim
+    real(wp_), dimension(:,:), allocatable, intent(out) :: psin
+    integer, optional, intent(in) :: ipsinorm,idesc,ifreefmt,unit
+! local variables
+    integer, parameter :: indef=0,iddef=1,iffdef=0
+    integer :: in,id,iffmt,u,idum,i,j,nr,nz,nbnd,nlim
+    character(len=48) :: string
+    real(wp_) :: dr,dz,dps,rleft,zmid,zleft,xdum,psiedge,psiaxis
+    
+! set default values if optional arguments are absent
+    in=indef; id=iddef; iffmt=iffdef
+    if(present(ipsinorm)) in=ipsinorm
+    if(present(idesc)) id=idesc
+    if(present(ifreefmt)) iffmt=ifreefmt
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+
+! open G EQDSK file (see http://fusion.gat.com/efit/g_eqdsk.html)
+    open(file=trim(filenm),status='old',unit=u)
+
+! get size of main arrays and allocate them
+    if (id==1) then
+      read (u,'(a48,3i4)') string,idum,nr,nz
+    else
+      read (u,*) nr,nz
+    end if
+    if (allocated(rv)) deallocate(rv)
+    if (allocated(zv)) deallocate(zv)
+    if (allocated(psin)) deallocate(psin)
+    if (allocated(psinr)) deallocate(psinr)
+    if (allocated(fpol)) deallocate(fpol)
+    if (allocated(q)) deallocate(q)
+    allocate(rv(nr),zv(nz),psin(nr,nz),psinr(nr),fpol(nr),q(nr))
+
+! store 0D data and main arrays
+    if (iffmt==1) then
+      read (u,*) dr,dz,rvac,rleft,zmid
+      read (u,*) rax,zax,psiaxis,psiedge,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) xdum,xdum,xdum,xdum,xdum
+      read (u,*) (fpol(i),i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) (xdum,i=1,nr)
+      read (u,*) ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,*) (q(i),i=1,nr)
+    else
+      read (u,'(5e16.9)') dr,dz,rvac,rleft,zmid
+      read (u,'(5e16.9)') rax,zax,psiaxis,psiedge,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') xdum,xdum,xdum,xdum,xdum
+      read (u,'(5e16.9)') (fpol(i),i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') (xdum,i=1,nr)
+      read (u,'(5e16.9)') ((psin(i,j),i=1,nr),j=1,nz)
+      read (u,'(5e16.9)') (q(i),i=1,nr)
+    end if
+
+! get size of boundary and limiter arrays and allocate them
+    read (u,*) nbnd,nlim
+    if (allocated(rbnd)) deallocate(rbnd)
+    if (allocated(zbnd)) deallocate(zbnd)
+    if (allocated(rlim)) deallocate(rlim)
+    if (allocated(zlim)) deallocate(zlim)
+
+! store boundary and limiter data
+    if(nbnd>0) then
+      allocate(rbnd(nbnd),zbnd(nbnd))
+      if (iffmt==1) then
+        read(u,*) (rbnd(i),zbnd(i),i=1,nbnd)
+      else
+        read(u,'(5e16.9)') (rbnd(i),zbnd(i),i=1,nbnd)
+      end if
+    end if
+    if(nlim>0) then
+      allocate(rlim(nlim),zlim(nlim))
+      if (iffmt==1) then
+        read(u,*) (rlim(i),zlim(i),i=1,nlim)
+      else
+        read(u,'(5e16.9)') (rlim(i),zlim(i),i=1,nlim)
+      end if
+    end if
+
+! reading of G EQDSK file completed    
+    close(u)
+    
+! build rv,zv,psinr arrays and normalize psin
+    zleft=zmid-0.5_wp_*dz
+    dr=dr/(nr-1)
+    dz=dz/(nz-1)
+    dps=one/(nr-1)
+    do i=1,nr
+      psinr(i)=(i-1)*dps
+      rv(i)=rleft+(i-1)*dr
+    end do
+    do i=1,nz
+      zv(i)=zleft+(i-1)*dz
+    end do
+    psia=psiedge-psiaxis
+    if(in==0) psin=(psin-psiaxis)/psia
+  end subroutine read_eqdsk
+
+  subroutine read_equil_an(filenm,rv,zv,fpol,q,unit)
+    use utils, only : get_free_unit
+    implicit none
+! arguments
+    character(len=*), intent(in) :: filenm
+    integer, optional, intent(in) :: unit
+!    integer, intent(in) :: isgnbphi
+!    real(wp_), intent(in) :: factb
+!    real(wp_), intent(out) :: rvac,rax,zax
+    real(wp_), dimension(:), allocatable, intent(out) :: rv,zv,fpol,q
+! local variables
+    integer :: u
+    real(wp_) :: rr0m,zr0m,rpam,b0,q0,qa,alq !,rcen,btrcen
+
+    if (present(unit)) then
+      u=unit
+    else
+      u=get_free_unit()
+    end if
+    open(file=trim(filenm),status='old',unit=u)
+    read(u,*) rr0m,zr0m,rpam
+    read(u,*) b0
+    read(u,*) q0,qa,alq
+!    rcen=rr0m
+!    btrcen=b0
+!    b0=isgnbphi*b0*factb
+!    rvac=rr0m
+!    rax=rr0m
+!    zax=zr0m
+!    zbmin=(zr0-rpa)/1.0e2_wp_
+!    zbmax=(zr0+rpa)/1.0e2_wp_
+    if(allocated(rv)) deallocate(rv)
+    if(allocated(zv)) deallocate(zv)
+    if(allocated(fpol)) deallocate(fpol)
+    if(allocated(q)) deallocate(q)
+    allocate(rv(2),zv(1),fpol(1),q(3))
+    rv(1)=rr0m
+    rv(2)=rpam
+    zv(1)=zr0m
+    fpol(1)=b0*rr0m
+    q(1)=q0
+    q(2)=qa
+    q(3)=alq
+    close(u)
+  end subroutine read_equil_an
+  
+  subroutine change_cocos(psia,fpol,q,cocosin,cocosout,ierr)
+    use const_and_precisions, only : zero,one,pi
+    implicit none
+! arguments
+    real(wp_), intent(inout) :: psia
+    real(wp_), dimension(:), intent(inout) :: fpol,q
+    integer, intent(in) :: cocosin, cocosout
+!    real(wp_), intent(out) :: isign,bsign
+    integer, intent(out), optional :: ierr
+! local variables
+    real(wp_) :: isign,bsign
+    integer :: exp2pi,exp2piout
+    logical :: phiccw,psiincr,qpos,phiccwout,psiincrout,qposout
+
+    call decode_cocos(cocosin,exp2pi,phiccw,psiincr,qpos)
+    call decode_cocos(cocosout,exp2piout,phiccwout,psiincrout,qposout)
+
+!   check sign consistency
+    isign=sign(one,psia)
+    if (.not.psiincr) isign=-isign
+    bsign=sign(one,fpol(size(fpol)))
+    if (qpos.neqv.isign*bsign*q(size(q))>zero) then
+!     warning: sign inconsistency found among q, Ipla and Bref
+      q=-q
+      if(present(ierr)) ierr=1
+    else
+      if(present(ierr)) ierr=0
+    end if
+
+! convert cocosin to cocosout
+    if (phiccw.neqv.phiccwout) then
+!     opposite direction of toroidal angle phi in cocosin and cocosout
+!      bsign=-bsign
+!      isign=-isign
+      fpol=-fpol
+    end if
+!   q has opposite sign for given sign of Bphi*Ip
+    if (qpos .neqv. qposout) q=-q
+!   psi and Ip signs don't change accordingly
+    if ((phiccw.eqv.phiccwout) .neqv. (psiincr.eqv.psiincrout)) psia=-psia
+!   convert Wb to Wb/rad or viceversa
+    psia=psia*(2.0_wp_*pi)**(exp2piout-exp2pi)
+  end subroutine change_cocos
+
+  subroutine decode_cocos(cocos,exp2pi,phiccw,psiincr,qpos)
+    implicit none
+    integer, intent(in) :: cocos
+    integer, intent(out) :: exp2pi
+    logical, intent(out) :: phiccw,psiincr,qpos
+    integer :: cmod10,cmod4
+
+    cmod10=mod(cocos,10)
+    cmod4=mod(cmod10,4)
+!  cocos>10 psi in Wb, cocos<10 psi in Wb/rad
+    exp2pi=(cocos-cmod10)/10
+!  cocos mod 10 = 1,3,5,7: toroidal angle phi increasing CCW
+    phiccw=(mod(cmod10,2)==1)
+!  cocos mod 10 = 1,2,5,6: psi increasing with positive Ip
+    psiincr=(cmod4==1 .or. cmod4==2)
+!  cocos mod 10 = 1,2,7,8: q positive for positive Bphi*Ip
+    qpos=(cmod10<3 .or. cmod10>6)
+  end subroutine decode_cocos
+
+  subroutine eq_scal(psia,fpol,isign,bsign,factb)
+    use const_and_precisions, only : one
+    implicit none
+    real(wp_), intent(inout) :: psia
+    integer, intent(inout) :: isign,bsign
+    real(wp_), dimension(:), intent(inout) :: fpol
+    real(wp_), intent(in) :: factb
+    
+    ! isign and bsign ignored on input if equal to 0
+    ! used to assign the direction of Bphi and Ipla BEFORE scaling
+    ! cocos=3 assumed: CCW direction is >0
+    ! Bphi and Ipla scaled by the same factor factb to keep q unchanged
+    ! factb<0 reverses the directions of Bphi and Ipla
+    if(isign/=0) psia=sign(psia,real(-isign,wp_))
+    if(bsign/=0 .and. bsign*fpol(size(fpol))<0) fpol=-fpol
+    psia=psia*factb
+    fpol=fpol*factb
+    isign=int(sign(one,-psia))
+    bsign=int(sign(one,fpol(size(fpol))))
+  end subroutine eq_scal
+
+  subroutine set_eqspl(rv,zv,psin,psiwbrad,psinr,fpol,sspl,ssfp, &
+      r0,rax,zax,rbnd,zbnd,ixp)
+    use const_and_precisions, only : zero,one
+    use dierckx, only : regrid,coeff_parder,curfit,splev
+    use utils, only : vmaxmin,vmaxmini
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(in) :: rv,zv,psinr,fpol
+    real(wp_), dimension(:,:), intent(in) :: psin
+    real(wp_), intent(in) :: psiwbrad
+    real(wp_), intent(inout) :: sspl,ssfp
+    real(wp_), intent(in), optional :: r0,rax,zax
+    real(wp_), dimension(:), intent(in), optional :: rbnd,zbnd
+    integer, intent(in), optional :: ixp
+! local constants
+    integer, parameter :: iopt=0
+! local variables
+    integer :: liwrk,lwrk,lw10,lw01,lw20,lw02,lw11,lwrkf
+    integer :: nr,nz,nrest,nzest,npsest,nrz,npsi,nbnd,ibinf,ibsup
+    real(wp_) :: fp,rax0,zax0,psinoptmp,psinxptmp,rbmin,rbmax,rbinf,rbsup,r1,z1
+    real(wp_), dimension(1) :: fpoli
+    real(wp_), dimension(:), allocatable :: fvpsi,wf,wrk
+    integer, dimension(:), allocatable :: iwrk
+    integer :: ier,ixploc,info
+!
+!   compute array sizes and prepare working space arrays
+    nr=size(rv)
+    nz=size(zv)
+    nrz=nr*nz
+    nrest=nr+ksplp
+    nzest=nz+ksplp
+    lwrk=4+nrest*nz+(nrest+nzest)*(2*kspl+5)+(nr+nz)*ksplp+max(nz,nrest)
+    liwrk=nz+nr+nrest+nzest+kspl
+!
+    npsi=size(psinr)
+    npsest=npsi+ksplp
+    lwrkf=npsi*ksplp+npsest*(7+3*kspl)
+!
+    allocate(wrk(max(lwrk,lwrkf)),iwrk(max(liwrk,npsest)))
+!
+!   spline fitting/interpolation of psin(i,j) and derivatives
+!
+!   length in m !!!
+!
+    rmnm=rv(1)
+    rmxm=rv(nr)
+    zmnm=zv(1)
+    zmxm=zv(nz)
+!   allocate knots and spline coefficients arrays
+    if (allocated(tr)) deallocate(tr)
+    if (allocated(tz)) deallocate(tz)
+    allocate(tr(nrest),tz(nzest),cceq(nrz))
+!   allocate work arrays
+!   reshape 2D psi array to 1D (transposed) array and compute spline coeffs
+    allocate(fvpsi(nrz))
+    fvpsi=reshape(transpose(psin),(/nrz/))
+    call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,         &
+                kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp,   &
+                wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
+!   if ier=-1 data are re-fitted using sspl=0    
+    if(ier==-1) then
+      sspl=0.0_wp_
+      call regrid(iopt,nr,rv,nz,zv,fvpsi,rmnm,rmxm,zmnm,zmxm,       &
+                  kspl,kspl,sspl,nrest,nzest,nsr,tr,nsz,tz,cceq,fp, &
+                  wrk(1:lwrk),lwrk,iwrk(1:liwrk),liwrk,ier)
+    end if
+    deallocate(fvpsi)
+!   compute spline coefficients for psi partial derivatives
+    lw10 = nr*(ksplp-1) + nz*ksplp     + nrz
+    lw01 = nr*ksplp     + nz*(ksplp-1) + nrz
+    lw20 = nr*(ksplp-2) + nz*ksplp     + nrz
+    lw02 = nr*ksplp     + nz*(ksplp-2) + nrz
+    lw11 = nr*(ksplp-1) + nz*(ksplp-1) + nrz
+    if (allocated(cceq10)) deallocate(cceq10)
+    if (allocated(cceq01)) deallocate(cceq01)
+    if (allocated(cceq20)) deallocate(cceq20)
+    if (allocated(cceq02)) deallocate(cceq02)
+    if (allocated(cceq11)) deallocate(cceq11)
+    allocate(cceq10(lw10),cceq01(lw01),cceq20(lw20),cceq02(lw02),cceq11(lw11))
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,0,cceq10,lw10,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,1,cceq01,lw01,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,2,0,cceq20,lw20,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,0,2,cceq02,lw02,ier)
+    call coeff_parder(tr,nsr,tz,nsz,cceq,kspl,kspl,1,1,cceq11,lw11,ier)
+!
+!   spline interpolation of fpol(i)
+!
+!   allocate knots and spline coefficients arrays
+    if (allocated(tfp)) deallocate(tfp)
+    if (allocated(cfp)) deallocate(cfp)
+    allocate(tfp(npsest),cfp(npsest))
+    allocate(wf(npsi))
+    wf(1:npsi-1)=one
+    wf(npsi)=1.0e2_wp_
+    call curfit(iopt,npsi,psinr,fpol,wf,zero,one,kspl,ssfp,npsest,nsf, &
+                tfp,cfp,fp,wrk(1:lwrkf),lwrkf,iwrk(1:npsest),ier)
+    call splev(tfp,nsf,cfp,kspl,psinr(npsi:npsi),fpoli,1,ier)
+!   set vacuum value used outside 0<=psin<=1 range
+    fpolas=fpoli(1)
+!   free temporary arrays
+    deallocate(wrk,iwrk,wf)
+!
+!   re-normalize psi after spline computation
+!
+!   start with un-corrected psi
+!
+    psia=psiwbrad
+    psinop=0.0_wp_
+    psiant=1.0_wp_
+!
+!   use provided boundary to set an initial guess for the search of O/X points
+!
+    nbnd=0
+    if(present(rbnd).and.present(zbnd)) then
+      nbnd=min(size(rbnd),size(zbnd))
+    end if
+    if (nbnd>0) then
+      call vmaxmini(zbnd,nbnd,zbinf,zbsup,ibinf,ibsup)
+      rbinf=rbnd(ibinf)
+      rbsup=rbnd(ibsup)
+      call vmaxmin(rbnd,nbnd,rbmin,rbmax)
+    else
+      zbinf=zv(2)
+      zbsup=zv(nz-1)
+      rbinf=rv((nr+1)/2)
+      rbsup=rbinf
+      rbmin=rv(2)
+      rbmax=rv(nr-1)
+    end if
+!
+!   search for exact location of the magnetic axis
+!
+    if(present(rax)) then
+      rax0=rax
+    else
+      rax0=0.5_wp_*(rbmin+rbmax)
+    end if
+    if(present(zax)) then
+      zax0=zax
+    else
+      zax0=0.5_wp_*(zbinf+zbsup)
+    end if
+    call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
+    print'(a,2f8.4,es12.5)','O-point',rmaxis,zmaxis,psinoptmp
+!
+!   search for X-point   if ixp not = 0
+!
+    if(present(ixp)) then
+      ixploc=ixp
+    else
+      ixploc=0
+    end if
+    if(ixploc/=0) then
+      if(ixploc<0) then
+        call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
+        if(psinxptmp/=-1.0_wp_) then
+          print'(a,2f8.4,es12.5)','X-point',r1,z1,psinxptmp
+          zbinf=z1
+          psinop=psinoptmp
+          psiant=psinxptmp-psinop
+          call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
+          zbsup=z1
+        else
+          ixploc=0
+        end if
+      else
+        call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
+        if(psinxptmp.ne.-1.0_wp_) then
+          print'(a,2f8.4,e16.8)','X-point',r1,z1,psinxptmp
+          zbsup=z1
+          psinop=psinoptmp
+          psiant=psinxptmp-psinop
+          call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
+          zbinf=z1
+        else
+          ixploc=0
+        end if
+      end if
+    end if
+
+    if (ixploc==0) then
+      psinop=psinoptmp
+      psiant=one-psinop
+!   find upper horizontal tangent point
+      call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
+      zbsup=z1
+      rbsup=r1
+!   find lower horizontal tangent point
+      call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
+      zbinf=z1
+      rbinf=r1
+      print'(a,4f8.4)','no X-point  ',rbinf,zbinf,rbsup,zbsup
+    end if
+    print*,' '
+!
+!   save Bt value on axis (required in flux_average and used in Jcd def)
+!   and vacuum value B0 at ref. radius R0 (used in Jcd_astra def)
+!
+    call equinum_fpol(0.0_wp_,btaxis)
+    btaxis=btaxis/rmaxis
+    if(present(r0)) then
+      btrcen=fpolas/r0
+      rcen=r0
+    else
+      btrcen=fpolas/rmaxis
+      rcen=rmaxis
+    end if
+    print'(a,f8.4)','BT_centr= ',btrcen
+    print'(a,f8.4)','BT_axis = ',btaxis
+  end subroutine set_eqspl
+
+  subroutine unset_eqspl
+    implicit none
+    if(allocated(tr)) deallocate(tr)
+    if(allocated(tz)) deallocate(tz)
+    if(allocated(tfp)) deallocate(tfp)
+    if(allocated(cfp)) deallocate(cfp)
+    if(allocated(cceq)) deallocate(cceq)
+    if(allocated(cceq01)) deallocate(cceq01)
+    if(allocated(cceq10)) deallocate(cceq10)
+    if(allocated(cceq02)) deallocate(cceq02)
+    if(allocated(cceq20)) deallocate(cceq20)
+    if(allocated(cceq11)) deallocate(cceq11)
+    nsr=0
+    nsz=0
+    nsf=0
+  end subroutine unset_eqspl
+
+  subroutine equinum_psi(rpsim,zpsim,psinv,dpsidr,dpsidz, &
+        ddpsidrr,ddpsidzz,ddpsidrz)
+    use dierckx, only : fpbisp
+    implicit none
+! local constants
+    integer, parameter :: lwrk=8,liwrk=2
+! arguments
+    real(wp_), intent(in) :: rpsim,zpsim
+    real(wp_), intent(out), optional :: psinv,dpsidr,dpsidz, &
+        ddpsidrr,ddpsidzz,ddpsidrz
+
+! local variables
+    integer, dimension(liwrk) :: iwrk
+    real(wp_), dimension(1) :: rrs,zzs,ffspl
+    real(wp_), dimension(lwrk) :: wrk
+!
+!     here lengths are measured in meters
+!
+    if (rpsim.le.rmxm .and. rpsim.ge.rmnm .and. &
+        zpsim.le.zmxm .and. zpsim.ge.zmnm) then
+
+      if (present(psinv)) then
+        rrs(1)=rpsim
+        zzs(1)=zpsim
+        call fpbisp(tr,nsr,tz,nsz,cceq,3,3,rrs,1,zzs,1,ffspl, &
+                   wrk(1),wrk(5),iwrk(1),iwrk(2))
+        psinv=(ffspl(1)-psinop)/psiant
+      end if
+      if (present(dpsidr)) then
+        call sub_derpsi(rpsim,zpsim,1,0,dpsidr,cceq10)
+      end if
+      if (present(dpsidz)) then
+        call sub_derpsi(rpsim,zpsim,0,1,dpsidz,cceq01)
+      end if
+      if (present(ddpsidrr)) then
+        call sub_derpsi(rpsim,zpsim,2,0,ddpsidrr,cceq20)
+      end if
+      if (present(ddpsidzz)) then
+        call sub_derpsi(rpsim,zpsim,0,2,ddpsidzz,cceq02)
+      end if
+      if (present(ddpsidrz)) then
+        call sub_derpsi(rpsim,zpsim,1,1,ddpsidrz,cceq11)
+      end if
+    else
+      if(present(psinv)) psinv=-1.0_wp_
+      if(present(dpsidr)) dpsidr=0.0_wp_
+      if(present(dpsidz)) dpsidz=0.0_wp_
+      if(present(ddpsidrr)) ddpsidrr=0.0_wp_
+      if(present(ddpsidzz)) ddpsidzz=0.0_wp_
+      if(present(ddpsidrz)) ddpsidrz=0.0_wp_
+    end if
+  end subroutine equinum_psi
+
+  subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc)
+    use dierckx, only : fpbisp
+    implicit none
+! arguments
+    real(wp_), intent(in) :: rpsim,zpsim
+    integer, intent(in) :: nur,nuz
+    real(wp_), intent(out) :: derpsi
+    real(wp_), dimension(:) :: cc
+! local variables
+    integer, dimension(1) :: iwrkr,iwrkz
+    real(wp_), dimension(1) :: rrs,zzs,ffspl
+    real(wp_), dimension(1,ksplp) :: wrkr
+    real(wp_), dimension(1,ksplp) :: wrkz
+
+    rrs(1)=rpsim
+    zzs(1)=zpsim
+    call fpbisp(tr(nur+1),nsr-2*nur,tz(nuz+1),nsz-2*nuz,cc,kspl-nur,kspl-nuz, &
+         rrs,1,zzs,1,ffspl,wrkr,wrkz,iwrkr,iwrkz)
+    derpsi=ffspl(1)*psia
+  end subroutine sub_derpsi
+
+  subroutine equinum_fpol(psinv,fpolv,dfpv)
+    use dierckx, only : splev,splder
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psinv
+    real(wp_), intent(out) :: fpolv
+    real(wp_), intent(out), optional :: dfpv
+! local variables
+    integer :: ier
+    real(wp_), dimension(1) :: rrs,ffspl
+    real(wp_), dimension(nsf) :: wrkfd
+!      
+    if(psinv.le.1.0_wp_.and.psinv.ge.0.0_wp_) then
+      rrs(1)=psinv
+      call splev(tfp,nsf,cfp,3,rrs,ffspl,1,ier)
+      fpolv=ffspl(1)
+      if(present(dfpv)) then       
+        call splder(tfp,nsf,cfp,3,1,rrs,ffspl,1,wrkfd,ier)
+        dfpv=ffspl(1)/psia
+      end if
+    else
+      fpolv=fpolas
+      if (present(dfpv)) dfpv=0._wp_
+    end if
+  end subroutine equinum_fpol
+
+  subroutine bfield(rpsim,zpsim,bphi,br,bz)
+    implicit none
+! arguments
+    real(wp_), intent(in) :: rpsim,zpsim
+    real(wp_), intent(out), optional :: bphi,br,bz
+! local variables
+    real(wp_) :: psin,fpol
+
+    call equinum_psi(rpsim,zpsim,psinv=bphi,dpsidr=bz,dpsidz=br)
+    if (present(bphi)) then
+      psin=bphi
+      call equinum_fpol(psin,fpol)
+      bphi=fpol/rpsim
+    end if
+    if (present(br)) br=-br/rpsim
+    if (present(bz)) bz= bz/rpsim
+  end subroutine bfield
+
+  subroutine setqphi_num(psinq,q,psia,rhotn)
+    use const_and_precisions, only : pi
+    use simplespline, only : difcs
+    implicit none
+!       arguments
+    real(wp_), dimension(:), intent(in) :: psinq,q
+    real(wp_), intent(in) :: psia
+    real(wp_), dimension(:), intent(out), optional :: rhotn
+!       local variables
+    real(wp_), dimension(size(q)) :: phit
+    real(wp_) :: dx
+    integer, parameter :: iopt=0
+    integer :: k,ier
+
+    nq=size(q)
+    if(allocated(psinr)) deallocate(psinr)
+    if(allocated(cq)) deallocate(cq)
+    allocate(psinr(nq),cq(nq,4))
+
+    psinr=psinq
+    call difcs(psinr,q,nq,iopt,cq,ier)
+!
+!   Toroidal flux phi = 2*pi*Integral q dpsi
+!
+    phit(1)=0.0_wp_
+    do k=1,nq-1
+      dx=psinr(k+1)-psinr(k)
+      phit(k+1)=phit(k) + dx*(cq(k,1)         + dx*(cq(k,2)/2.0_wp_ +  &
+                          dx*(cq(k,3)/3.0_wp_ + dx* cq(k,4)/4.0_wp_) ) )
+    end do
+    phitedge=phit(nq)
+    if(present(rhotn)) rhotn(1:nq)=sqrt(phit/phitedge)
+    phitedge=2*pi*psia*phitedge
+  end subroutine setqphi_num
+
+  subroutine unset_q
+    implicit none
+    
+    if(allocated(psinr)) deallocate(psinr)
+    if(allocated(cq)) deallocate(cq)
+    nq=0
+  end subroutine unset_q
+
+  function fq(psin)
+    use const_and_precisions, only : wp_
+    use simplespline, only :spli
+    use utils, only : locate
+    implicit none
+! arguments
+    real(wp_), intent(in) :: psin
+    real(wp_) :: fq
+! local variables
+    integer :: i
+    real(wp_) :: dps
+
+    call locate(psinr,nq,psin,i)
+    i=min(max(1,i),nq-1)
+    dps=psin-psinr(i)
+    fq=spli(cq,nq,i,dps)
+  end function fq
+
+  subroutine set_rhospl(rhop,rhot)
+    use simplespline, only : difcs
+    implicit none
+! arguments
+    real(wp_), dimension(:), intent(in) :: rhop, rhot
+! local variables
+    integer, parameter :: iopt=0
+    integer :: ier
+    
+    nrho=size(rhop)
+
+    if(allocated(rhopr)) deallocate(rhopr)
+    if(allocated(rhotr)) deallocate(rhotr)
+    if(allocated(crhop)) deallocate(crhop)
+    if(allocated(crhot)) deallocate(crhot)
+    allocate(rhopr(nrho),rhotr(nrho),crhop(nrho,4),crhot(nrho,4))
+    
+    rhopr=rhop
+    rhotr=rhot
+    call difcs(rhotr,rhopr,nrho,iopt,crhop,ier)
+    call difcs(rhopr,rhotr,nrho,iopt,crhot,ier)
+  end subroutine set_rhospl
+
+  subroutine unset_rhospl
+    implicit none
+    
+    if(allocated(rhopr)) deallocate(rhopr)
+    if(allocated(rhotr)) deallocate(rhotr)
+    if(allocated(crhop)) deallocate(crhop)
+    if(allocated(crhot)) deallocate(crhot)
+    nrho=0
+  end subroutine unset_rhospl
+
+  function frhopol(rhot)
+    use utils, only : locate
+    use simplespline, only : spli
+    implicit none
+!       arguments
+    real(wp_), intent(in) :: rhot
+    real(wp_) :: frhopol
+!       local variables
+    integer :: i
+    real(wp_) :: dr
+
+    call locate(rhotr,nrho,rhot,i)
+    i=min(max(1,i),nrho-1)
+    dr=rhot-rhotr(i)
+    frhopol=spli(crhop,nrho,i,dr)
+  end function frhopol
+
+  function frhotor(rhop)
+    use utils, only : locate
+    use simplespline, only : spli
+    implicit none
+!       arguments
+    real(wp_), intent(in) :: rhop
+    real(wp_) :: frhotor
+!       local variables
+    integer :: i
+    real(wp_) :: dr
+
+    call locate(rhopr,nrho,rhop,i)
+    i=min(max(1,i),nrho-1)
+    dr=rhop-rhopr(i)
+    frhotor=spli(crhot,nrho,i,dr)
+  end function frhotor
+
+  subroutine points_ox(rz,zz,rf,zf,psinvf,info)
+    use const_and_precisions, only : comp_eps
+    use minpack, only : hybrj1
+    implicit none
+! local constants
+    integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
+! arguments
+    real(wp_), intent(in) :: rz,zz
+    real(wp_), intent(out) :: rf,zf,psinvf
+    integer, intent(out) :: info
+! local variables
+    real(wp_) :: tol
+    real(wp_), dimension(n) :: xvec,fvec
+    real(wp_), dimension(lwa) :: wa
+    real(wp_), dimension(ldfjac,n) :: fjac
+
+    xvec(1)=rz
+    xvec(2)=zz
+    tol = sqrt(comp_eps)
+    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, &
+          '  O/X coord.',xvec
+      end if
+      rf=xvec(1)
+      zf=xvec(2)
+      call equinum_psi(rf,zf,psinvf)
+  end subroutine points_ox
+
+  subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
+    implicit none
+! arguments
+    integer, intent(in) :: n,iflag,ldfjac
+    real(wp_), dimension(n), intent(in) :: x
+    real(wp_), dimension(n), intent(inout) :: fvec
+    real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
+! local variables
+    real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
+!
+    select case(iflag)
+    case(1)
+      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
+      fvec(1) = dpsidr/psia
+      fvec(2) = dpsidz/psia    
+    case(2)
+      call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, &
+        ddpsidrz=ddpsidrz)
+      fjac(1,1) = ddpsidrr/psia
+      fjac(1,2) = ddpsidrz/psia
+      fjac(2,1) = ddpsidrz/psia
+      fjac(2,2) = ddpsidzz/psia
+    case default
+      print*,'iflag undefined'
+    end select
+  end subroutine fcnox
+
+  subroutine points_tgo(rz,zz,rf,zf,psin0,info)
+    use const_and_precisions, only : comp_eps
+    use minpack, only : hybrj1mv
+    implicit none
+! local constants
+    integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
+! arguments
+    real(wp_), intent(in) :: rz,zz,psin0
+    real(wp_), intent(out) :: rf,zf
+    integer, intent(out) :: info
+! local variables
+    real(wp_) :: tol
+    real(wp_), dimension(n) :: xvec,fvec,f0
+    real(wp_), dimension(lwa) :: wa
+    real(wp_), dimension(ldfjac,n) :: fjac
+! common/external functions/variables
+
+    xvec(1)=rz
+    xvec(2)=zz
+    f0(1)=psin0
+    f0(2)=0.0_wp_
+    tol = sqrt(comp_eps)
+    call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
+    if(info.gt.1)  then
+      print'(a,i2,a,2f8.4)','    info subr points_tgo =',info, &
+        '  R,z coord.',xvec
+    end if
+    rf=xvec(1)
+    zf=xvec(2)
+  end
+
+  subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
+    use const_and_precisions, only : wp_
+    implicit none
+! arguments
+    integer, intent(in) :: n,ldfjac,iflag
+    real(wp_), dimension(n), intent(in) :: x,f0
+    real(wp_), dimension(n), intent(inout) :: fvec
+    real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
+! internal variables
+    real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
+
+    select case(iflag)
+    case(1)
+      call equinum_psi(x(1),x(2),psinv,dpsidr)
+      fvec(1) = psinv-f0(1)
+      fvec(2) = dpsidr/psia-f0(2)
+    case(2)
+      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, &
+        ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
+      fjac(1,1) = dpsidr/psia
+      fjac(1,2) = dpsidz/psia
+      fjac(2,1) = ddpsidrr/psia
+      fjac(2,2) = ddpsidrz/psia
+    case default
+      print*,'iflag undefined'
+    end select
+  end subroutine fcntgo
+
+  subroutine set_equian(rax,zax,a,bax,qax,q1,qexp,n)
+    use const_and_precisions, only : pi,zero,one
+    implicit none
+! arguments
+    real(wp_), intent(in) :: rax,zax,a,bax,qax,q1,qexp
+    integer, intent(in), optional :: n
+! local variables
+    integer, parameter :: nqdef=101
+    integer :: i
+    real(wp_) :: dr,fq0,fq1,qq,res,rn
+    real(wp_), dimension(:), allocatable :: rhotn,rhopn
+
+    btaxis=bax
+    rmaxis=rax
+    zmaxis=zax
+    btrcen=bax
+    rcen=rax
+    aminor=a
+    zbinf=zmaxis-a
+    zbsup=zmaxis+a
+    q0=qax
+    qa=q1
+    alq=qexp
+
+    if (present(n)) then
+      nq=n
+    else
+      nq=nqdef
+    end if
+    if (allocated(psinr)) deallocate(psinr)
+    allocate(psinr(nq),rhotn(nq),rhopn(nq))
+    dr=one/(nq-1)
+    rhotn(1)=zero
+    psinr(1)=zero
+    res=zero
+    fq0=zero
+    do i=2,n
+      rn=(i-1)*dr
+      qq=q0+(q1-q0)*rn**qexp
+      fq1=rn/qq
+      res=res+0.5_wp_*(fq1+fq0)*dr
+      fq0=fq1
+      rhotn(i)=rn
+      psinr(i)=res
+    end do
+    phitedge=btaxis*aminor**2 ! temporary
+    psia=res*phitedge
+    phitedge=pi*phitedge    ! final
+    psinr=psinr/res
+    rhopn=sqrt(psinr)
+    call set_rhospl(rhopn,rhotn)
+  end subroutine set_equian
+
+  subroutine equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, &
+      ddpsidrr,ddpsidzz,ddpsidrz)
+    use const_and_precisions, only : wp_
+    implicit none
+! arguments
+    real(wp_), intent(in) :: rrm,zzm
+    real(wp_), intent(out), optional :: psinv,fpolv,dfpv,dpsidr,dpsidz, &
+      ddpsidrr,ddpsidzz,ddpsidrz
+! local variables
+    integer :: sgn
+    real(wp_) :: cst,dpsidrp,d2psidrp,dqq,qq,rn,rpm,snt,rhop,rhot
+!    real(wp_) :: frhopol
+
+!     simple model for equilibrium: large aspect ratio
+!     outside plasma: analytical continuation, not solution Maxwell eqs
+
+    rpm=sqrt((rrm-rmaxis)**2+(zzm-zmaxis)**2)  !!! rpm==rho_tor[m], rn=rho_tor_norm
+    rn=rpm/aminor
+
+    snt=0.0_wp_
+    cst=1.0_wp_
+    if (rpm > 0.0_wp_) then
+      snt=(zzm-zmaxis)/rpm
+      cst=(rrm-rmaxis)/rpm
+    end if
+
+    if (present(psinv)) then
+      rhot=rn
+      if(rn <= 1.0_wp_) then
+        rhop=frhopol(rhot)
+        psinv=rhop*rhop
+      else
+        psinv=1.0_wp_+btaxis*aminor**2/2.0_wp_/psia/qa*(rn*rn-1.0_wp_)
+        rhop=sqrt(psinv)
+      end if
+    end if
+
+    if(rn <= 1.0_wp_) then
+      qq=q0+(qa-q0)*rn**alq
+      dpsidrp=-btaxis*aminor*rn/qq
+      dqq=alq*(qa-q0)*rn**(alq-1.0_wp_)
+      d2psidrp=-btaxis*(1.0_wp_-rn*dqq/qq)/qq
+    else
+      dpsidrp=-btaxis*aminor/qa*rn
+      d2psidrp=-btaxis/qa
+    end if
+
+    if(present(fpolv)) fpolv=btaxis*rmaxis
+    if(present(dfpv))  dfpv=0.0_wp_      
+
+    if(present(dpsidr)) dpsidr=dpsidrp*cst
+    if(present(dpsidz)) dpsidz=dpsidrp*snt
+    if(present(ddpsidrr)) ddpsidrr=dpsidrp*snt**2/rpm+cst**2*d2psidrp
+    if(present(ddpsidrz)) ddpsidrz=cst*snt*(d2psidrp-dpsidrp/rpm)
+    if(present(ddpsidzz)) ddpsidzz=dpsidrp*cst**2/rpm+snt**2*d2psidrp
+
+  end subroutine equian
+
+end module equilibrium

src/graydata_anequil.f90

@@ -0,0 +1,12 @@
+module graydata_anequil
+  use const_and_precisions, only : wp_
+  implicit none
+  
+  real(wp_), save :: dens0,aln1,aln2
+  real(wp_), save :: te0,dte0,alt1,alt2
+  real(wp_), save :: rr0,zr0,rpa
+  real(wp_), save :: b0,rr0m,zr0m,rpam
+  real(wp_), save :: q0,qa,alq
+  real(wp_), save :: zeffan
+
+end module graydata_anequil

src/graydata_flags.f90

@@ -0,0 +1,14 @@
+module graydata_flags
+  use const_and_precisions, only : wp_
+  implicit none
+  
+  character*255, save :: filenmeqq,filenmprf,filenmbm
+  real(wp_), save :: sspl, dst
+  integer, save :: ibeam,iequil,ixp,iprof
+  integer, save :: iwarm,ilarm,imx,ieccd,ipec,idst
+  integer, save :: igrad,ipass
+  integer, save :: ipsinorm,iscal,icocos
+  integer, save :: nnd,istpr0,istpl0,ipol
+  integer, save :: neqdsk,nprof
+  
+end module graydata_flags

src/graydata_par.f90

@@ -0,0 +1,10 @@
+module graydata_par
+  use const_and_precisions, only : wp_
+  implicit none
+  
+  real(wp_), save :: rwmax,rwallm
+  real(wp_), save :: psipol0,chipol0
+  real(wp_), save :: factb,factt,factn
+  real(wp_), save :: sgnbphi,sgniphi
+
+end module graydata_par

src/green_func_p.f90

@@ -1,437 +0,0 @@
-!########################################################################
-
- MODULE green_func_p
-
-!########################################################################
-!
-! The module contains few subroutines which are requested to calculate 
-! the current drive value by adjoint approach
-!
-!########################################################################
- USE const_and_precisions
-!-------
- IMPLICIT NONE
- CHARACTER(Len=1), PRIVATE :: adj_appr(6)    ! adjoint approach switcher
-!-------
- REAL(wp_), PRIVATE :: r2,q2,gp1,Rfactor
-!-------
- REAL(wp_), PRIVATE, PARAMETER :: delta = 1e-4   ! border for recalculation
-!------- for N.M. subroutines (variational principle) -------
- REAL(wp_), PRIVATE :: sfd(1:4)
- INTEGER,   PRIVATE, PARAMETER :: nre = 2        ! order of rel. correct.
- REAL(wp_), PRIVATE, PARAMETER :: vp_mee(0:4,0:4,0:2) = &
-            RESHAPE((/0.0,      0.0,      0.0,      0.0,      0.0,      &
-                      0.0,      0.184875, 0.484304, 1.06069,  2.26175,  &
-                      0.0,      0.484304, 1.41421,  3.38514,  7.77817,  &
-                      0.0,      1.06069,  3.38514,  8.73232,  21.4005,  &
-                      0.0,      2.26175,  7.77817,  21.4005,  55.5079,  &
-                      !                                                 &
-                      0.0,     -1.33059,-2.57431, -5.07771, -10.3884,   &
-                     -0.846284,-1.46337, -1.4941, -0.799288, 2.57505,   &
-                     -1.1601,  -1.4941,   2.25114,  14.159,   50.0534,  &
-                     -1.69257, -0.799288, 14.159,   61.4168,  204.389,  &
-                     -2.61022,  2.57505,  50.0534,  204.389,  683.756,  &
-                      !                                                 &
-                      0.0,      2.62498,  0.985392,-5.57449, -27.683,   &
-                      0.0,      3.45785,  5.10096,  9.34463,  22.9831,  &
-                     -0.652555, 5.10096,  20.5135,  75.8022,  268.944,  &
-                     -2.11571,  9.34463,  75.8022,  330.42,   1248.69,  &
-                     -5.38358,  22.9831,  268.944,  1248.69,  4876.48/),&
-                    (/5,5,3/))
- REAL(wp_), PRIVATE, PARAMETER :: vp_mei(0:4,0:4,0:2) = &
-            RESHAPE((/0.0,     0.886227, 1.0,      1.32934,  2.0,       &
-                      0.886227,1.0,      1.32934,  2.0,      3.32335,   &
-                      1.0,     1.32934,  2.0,      3.32335,  6.0,       &
-                      1.32934, 2.0,      3.32335,  6.0,      11.6317,   &
-                      2.0,     3.32335,  6.0,      11.6317,  24.0,      &
-                      !                                                 &
-                      0.0,     0.332335, 1.0,      2.49251,  6.0,       &
-                      1.66168, 1.0,      2.49251,  6.0,      14.5397,   &
-                      3.0,     2.49251,  6.0,      14.5397,  36.0,      &
-                      5.81586, 6.0,      14.5397,  36.0,     91.5999,   &
-                      12.0,    14.5397,  36.0,     91.5999,  240.0,     &
-                      !                                                 &
-                      0.0,    -0.103855, 0.0,      1.09047,  6.0,       &
-                      0.726983,0.0,      1.09047,  6.0,      24.5357,   &
-                      3.0,     1.09047,  6.0,      24.5357,  90.0,      &
-                      9.81427, 6.0,      24.5357,  90.0,     314.875,   &
-                      30.0,    24.5357,  90.0,     314.875,  1080.0 /), &
-                    (/5,5,3/))
- REAL(wp_), PRIVATE, PARAMETER :: vp_oee(0:4,0:4,0:2) = &
-            RESHAPE((/0.0,     0.56419,  0.707107, 1.0073,   1.59099,   &
-                      0.56419, 0.707107, 1.0073,   1.59099,  2.73981,   &
-                      0.707107,1.0073,   1.59099,  2.73981,  5.08233,   &
-                      1.0073,  1.59099,  2.73981,  5.08233,  10.0627,   &
-                      1.59099, 2.73981,  5.08233,  10.0627,  21.1138,   &
-                      !                                                 &
-                      0.0,     1.16832,  1.90035,  3.5758,   7.41357,   &
-                      2.17562, 1.90035,  3.5758,   7.41357,  16.4891,   &
-                      3.49134, 3.5758,   7.41357,  16.4891,  38.7611,   &
-                      6.31562, 7.41357,  16.4891,  38.7611,  95.4472,   &
-                      12.4959, 16.4891,  38.7611,  95.4472,  244.803,   &
-                      !                                                 &
-                      0.0,     2.65931,  4.64177,  9.6032,   22.6941,   &
-                      4.8652,  4.64177,  9.6032,   22.6941,  59.1437,   &
-                      9.51418, 9.6032,   22.6941,  59.1437,  165.282,   &
-                      21.061,  22.6941,  59.1437,  165.282,  485.785,   &
-                      50.8982, 59.1437,  165.282,  485.785,  1483.22/), &
-                    (/5,5,3/))
- REAL(wp_), PRIVATE, PARAMETER :: vp_g(0:4,0:2) = &
-            RESHAPE((/1.32934, 2.0,      3.32335,  6.0,      11.6317,   &
-                      2.49251, 0.0,      2.90793,  12.0,     39.2571,   &
-                      1.09047, 6.0,      11.45,    30.0,     98.9606/), &
-                    (/5,3/))
-!########################################################################
-
- CONTAINS
-
-!#######################################################################
- 
- SUBROUTINE Setup_SpitzFunc
-!=======================================================================
- IMPLICIT NONE
-!=======================================================================
- adj_appr(1) = 'l'         ! collisionless limit
-! adj_appr(1) = 'c'         ! collisional (classical) limit, w/o trap. part.
- adj_appr(2) = 'm'         ! momentum conservation
-! adj_appr(2) = 'h'         ! high-speed limit
-!--- 
- adj_appr(3) = 'l'         ! DO NOT CHANGE!
- adj_appr(4) = 'r'         ! DO NOT CHANGE!
- adj_appr(5) = 'v'         ! DO NOT CHANGE!
- adj_appr(6) = 'i'         ! DO NOT CHANGE!
-!=======================================================================
-!.....
-!=======================================================================
- RETURN
- END SUBROUTINE Setup_SpitzFunc
-
- 
- SUBROUTINE GenSpitzFunc(Te,Zeff,fc,u,q,gam, K,dKdu)
-
-!=======================================================================
-! Author:  N.B.Marushchenko
-! June 2005: as start point the subroutine of Ugo Gasparino (198?) 
-!            SpitzFunc() is taken and modified.
-!            1. adapted to the Fortran-95
-!            2. derivative of Spitzer function is added
-!            3. separation for 2 brunches is done: 
-!               1st is referenced as 'with conservation of the moment',
-!               2nd - as 'high speed limit'. 
-!               The last one is taken from the Lin-Liu formulation 
-!               (Phys.Plasmas 10 (2003) 4064) with K = F*fc.
-!               The asymptotical high speed limit (Taguchi-Fisch model)
-!               is also included as the reference case.
-! Feb. 2008: non-relativ. version is replaced by the relativistic one;
-!            the method is the the same, but the trial-function is 
-!            based on the relativistic formulation.
-!            The relativistic corrections for the collisional operator 
-!            up to the second order, i.e. (1/mu)**2, are applied.
-! Sep. 2008: generalized Spitzer function for arbitrary collisionality
-!            is implemented. The model is based on the concept of 
-!            the "effective trapped particles fraction". 
-!            The different.-integral kinetic equation for the generalized
-!            Spitzer function is produced with help of subroutines 
-!            ArbColl_TrappFract_Array and ArbColl_SpitzFunc_Array,
-!            where the subroutines of H. Maassberg are called).
-!========================================================================
-! Spitzer function with & w/o trapped particle effects is given by:
-!
-!    K(x) = x/gamma*(d1*x+d2*x^2+d4*x^3+d4*x^4),
-! 
-! where x = v/v_th and gamma=1 for non-relativistic version (Ugo),
-! or    x = p/p_th for relativistic version (N.M., February 2008).
-! Note, that somewhere the function F(x) instead of K(x) is applied,
-!
-!    F(x) = K(x)/fc.
-!
-! Numerical inversion of the 5x5 symmetric matrix obtained from the
-! generalized Spitzer problem (see paper of Taguchi for the equation
-! and paper of Hirshman for the variational approach bringing to the
-! matrix to be inverted).
-! 
-! The numerical method used is an improved elimination scheme
-! (Banachiewiczs-Cholesky-Crout method).
-! This method is particularly simple for symmetric matrix.
-! As a reference see "Mathematical Handbook" by Korn & Korn, p.635-636.
-! 
-! Refs.:  1. S.P. Hirshman, Phys. Fluids 23 (1980) 1238
-!         2. M. Rome' et al., Plasma Phys. Contr. Fus. 40 (1998) 511
-!         3. N.B. Marushchenko et al., Fusion Sci. Technol. 55 (2009) 180
-!========================================================================
-! INPUTS:
-!  u    - p/sqrt(2mT)
-!  q    - p/mc;
-!  gam  - relativistic factor;
-!  mu   - mc2/Te
-!  Zeff - effective charge;
-!  fc   - fraction of circulating particles.
-!
-! OUTPUTS:
-!   K   - Spitzer's function
-!  dKdu = dK/du, i.e. its derivative over normalized momentum
-!=======================================================================
- IMPLICIT NONE
- REAL(wp_), INTENT(in)  :: Te,Zeff,fc,u,q,gam
- REAL(wp_), INTENT(out) :: K,dKdu
- REAL(wp_) :: mu,gam1,gam2,gam3,w,dwdu
-!=======================================================================
-  K   = 0
- dKdu = 0
- IF (u < comp_eps) RETURN
-!---
- mu = mc2_/max(Te,1d-3)
- SELECT CASE(adj_appr(2))
-   CASE('m')  !--------------- momentum conservation ------------------!
-     gam1 = gam                                                        !
-     IF (adj_appr(4) == 'n') gam1 = 1                                  !
-     gam2 = gam1*gam1                                                  !
-     gam3 = gam1*gam2                                                  !
-      K   = u/gam1*u*(sfd(1)+u*(sfd(2)+u*(sfd(3)+u*sfd(4))))           !
-     dKdu = u/gam3*  (sfd(1)*(1+  gam2)+u*(sfd(2)*(1+2*gam2)+ &        !
-                   u*(sfd(3)*(1+3*gam2)+u* sfd(4)*(1+4*gam2))))        !
-   !--------------------- end momentum conservation -------------------!
-   CASE('h')  !---------------- high-speed-limit ----------------------!
-     IF (adj_appr(4) == 'n') THEN     !- non-relativ. asymptotic form -!
-       K   =   u**4 *fc/(Zeff+1+4*fc) !- (Taguchi-Fisch model)        -!
-      dKdu = 4*u**3 *fc/(Zeff+1+4*fc)                                  !
-     ELSEIF (adj_appr(4) == 'r') THEN !- relativistic, Lin-Liu form.  -!
-       CALL SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)          !
-     ENDIF                                                             !
-   CASE default   !----------------------------------------------------!
-     PRINT*,'GenSpitzFunc: WARNING! Spitzer function is not defined.'
-     RETURN
- END SELECT
-!=======================================================================
- RETURN
- END SUBROUTINE GenSpitzFunc
-
-!#######################################################################
-!#######################################################################
-!#######################################################################
-
- SUBROUTINE SpitzFuncCoeff(Te,Zeff,fc)
-!=======================================================================
-! Calculates the matrix coefficients required for the subroutine 
-! "GenSpitzFunc", where the Spitzer function is defined through the
-! variational principle.
-! 
-! Weakly relativistic (upgraded) version (10.09.2008). 
-! Apart of the non-relativistic matrix coefficients, taken from the
-! old subroutine of Ugo Gasparino, the relativistic correction written
-! as series in 1/mu^n (mu=mc2/T) powers is added. Two orders are taken 
-! into account, i.e. n=0,1,2.
-!
-! In this version, the coefficients "oee", i.e. Omega_ij, are formulated
-! for arbitrary collisionality. 
-!
-! INPUT VARIABLES:
-! rho  = sqrt(SS) with SS - flux-surface label (norm. magn. flux)
-! ne   - density, 1/m^3
-! Te   - temperature, keV
-! Zeff - effective charge
-! fc   - fraction of circulating particles
-!
-! OUTPUT VARIABLES (defined as a global ones): 
-! sfd(1),...,sfd(4) - coefficients of the polynomial expansion of the 
-!                 "Spitzer"-function (the same as in the Hirshman paper)
-!=======================================================================
- IMPLICIT NONE
- REAL(wp_), INTENT(in) :: Te,Zeff,fc
- INTEGER   :: n,i,j
- REAL(wp_) :: rtc,rtc1,mu,y,tn(1:nre)
- REAL(wp_) :: m(0:4,0:4),g(0:4)
- REAL(wp_) :: om(0:4,0:4)
- REAL(wp_) :: gam11,gam21,gam31,gam41,gam01, &
-                    gam22,gam32,gam42,gam02, &
-                          gam33,gam43,gam03, &
-                                gam44,gam04,gam00
- REAL(wp_) :: alp12,alp13,alp14,alp10, &
-                    alp23,alp24,alp20, &
-                          alp34,alp30,alp40
- REAL(wp_) :: bet0,bet1,bet2,bet3,bet4,d0
- LOGICAL   :: renew,rel,newTe,newne,newZ,newfc
- REAL(wp_), SAVE :: sfdx(1:4) = 0
- REAL(wp_), SAVE :: ne_old =-1, Te_old =-1, Zeff_old =-1, fc_old =-1
-!=======================================================================
- rel   = Te > 1
- newTe = abs(Te  -Te_old  ) > delta*Te
- newZ  = abs(Zeff-Zeff_old) > delta*Zeff
- newfc = abs(fc  -fc_old  ) > delta*fc
- SELECT CASE(adj_appr(1))
-   CASE ('l','c')
-     renew = (newTe .and. rel) .OR. newZ .OR. newfc
- END SELECT
-!---
- IF (.not.renew) THEN
-   sfd(:) = sfdx(:)
-   RETURN
- ENDIF
-!=======================================================================
- tn(:) = 0
- IF (adj_appr(4) == 'r') THEN
-   IF (nre > 0) THEN
-     mu = mc2_/max(Te,1d-3)
-     tn(1) = 1/mu
-     DO n=2,min(2,nre)
-       tn(n) = tn(n-1)/mu
-     ENDDO
-   ENDIF
- ENDIF
-!---
- SELECT CASE(adj_appr(1))
-   CASE ('l','c')     !---- both classical & collisionless limits  ----!
-     rtc = (1-fc)/fc;   rtc1 = rtc+1                                   !
-     !---                                                              !
-     DO i=0,4                                                          !
-       g(i) = vp_g(i,0)                                                !
-       DO n=1,min(2,nre)                                               !
-         g(i) = g(i) + tn(n)*vp_g(i,n)                                 !
-       ENDDO                                                           !
-       !---                                                            !
-       DO j=0,4                                                        !
-         IF (i == 0 .or. j == 0 .or. j >= i) THEN                      !
-             y =      vp_mee(i,j,0) + rtc *vp_oee(i,j,0) + &           !
-                                 Zeff*rtc1*vp_mei(i,j,0)               !
-           DO n=1,min(2,nre)                                           !
-             y = y + (vp_mee(i,j,n) + rtc *vp_oee(i,j,n) + &           !
-                                 Zeff*rtc1*vp_mei(i,j,n))*tn(n)        !
-           ENDDO                                                       !
-           m(i,j) = y                                                  !
-         ENDIF                                                         !
-       ENDDO                                                           !
-     ENDDO                                                             !
-     DO i=2,4                                                          !
-       DO j=1,i-1                                                      !
-         m(i,j) = m(j,i)                                               !
-       ENDDO                                                           !
-     ENDDO                                                             !
-     m(0,0) = 0                                                        !
-   CASE default       !------------------------------------------------!
-     PRINT*,'Green_Func: WARNING! Adjoint approach is not defined.'
-     RETURN
- END SELECT
-!=======================================================================
- gam11 = m(1,1)
- gam21 = m(2,1)
- gam31 = m(3,1)
- gam41 = m(4,1)
- gam01 = m(0,1)
-!
- alp12 = m(1,2)/m(1,1)
- alp13 = m(1,3)/m(1,1)
- alp14 = m(1,4)/m(1,1)
- alp10 = m(1,0)/m(1,1)
-! 
- gam22 = m(2,2)-gam21*alp12
- gam32 = m(3,2)-gam31*alp12
- gam42 = m(4,2)-gam41*alp12
- gam02 = m(0,2)-gam01*alp12
-!
- alp23 = gam32/gam22
- alp24 = gam42/gam22
- alp20 = gam02/gam22
-!
- gam33 = m(3,3)-gam31*alp13-gam32*alp23
- gam43 = m(4,3)-gam41*alp13-gam42*alp23
- gam03 = m(0,3)-gam01*alp13-gam02*alp23
-!
- alp34 = gam43/gam33
- alp30 = gam03/gam33
-!
- gam44 = m(4,4)-gam41*alp14-gam42*alp24-gam43*alp34
- gam04 = m(0,4)-gam01*alp14-gam02*alp24-gam03*alp34
-!
- alp40 = gam04/gam44
-!
- gam00 = m(0,0)-gam01*alp10-gam02*alp20-gam03*alp30-gam04*alp40
-!
- bet1 =  g(1)/m(1,1)
- bet2 = (g(2)-gam21*bet1)/gam22
- bet3 = (g(3)-gam31*bet1-gam32*bet2)/gam33
- bet4 = (g(4)-gam41*bet1-gam42*bet2-gam43*bet3)/gam44
- bet0 = (g(0)-gam01*bet1-gam02*bet2-gam03*bet3-gam04*bet4)/gam00
-!
- d0     = bet0
- sfd(4) = bet4-alp40*d0
- sfd(3) = bet3-alp30*d0-alp34*sfd(4)
- sfd(2) = bet2-alp20*d0-alp24*sfd(4)-alp23*sfd(3)
- sfd(1) = bet1-alp10*d0-alp14*sfd(4)-alp13*sfd(3)-alp12*sfd(2)
-!=======================================================================
- fc_old   = fc
- Te_old   = Te
- Zeff_old = Zeff
-!---
- sfdx(1:4) = sfd(1:4)
-!=======================================================================
- RETURN
- END SUBROUTINE SpitzFuncCoeff
-
-!#######################################################################
-!#######################################################################
-!#######################################################################
-
- SUBROUTINE SpitzFunc_HighSpeedLimit(Zeff,fc,u,q,gam, K,dKdu)
-!=======================================================================
-! Calculates the "Spitzer function" in high velocity limit, relativistic 
-! formulation: Lin-Liu et al., Phys.Pl. (2003),v10, 4064, Eq.(33).
-! 
-! Inputs:
-! Zeff - effective charge
-! fc   - fraction of circulating electrons
-! u    - p/(m*vte)
-! q    - p/mc
-! gam  - relativ. factor
-! 
-! Outputs:
-!  K   - Spitzer function
-! dKdu - its derivative
-!=======================================================================
- IMPLICIT NONE
- REAL(wp_), INTENT(in)  :: Zeff,fc,u,q,gam
- REAL(wp_), INTENT(out) :: K,dKdu
- INTEGER :: nfun
- REAL(8) :: gam2,err,flag,Integr
- REAL(8), PARAMETER :: a = 0d0, b = 1d0, rtol = 1d-4, atol = 1d-12
-!=======================================================================
- r2 = (1+Zeff)/fc   ! global parameter needed for integrand, HSL_f(t)
-!------------------
- IF (u < 1e-2) THEN
-    K   =   u**4/(r2+4)
-   dKdu = 4*u**3/(r2+4)
-   RETURN
- ENDIF
-!=======================================================================
- q2  = q*q       ! for the integrand, HSL_f
- gp1 = gam+1     ! ..
-!---
- CALL quanc8(HSL_f,zero,unit,atol,rtol,Integr,err,nfun,flag)
-!=======================================================================
- gam2 = gam*gam
-!---
-  K   = u**4 * Integr
- dKdu = (u/gam)**3 * (1-r2*gam2*Integr)
-!=======================================================================
- RETURN
- END SUBROUTINE SpitzFunc_HighSpeedLimit
-
-!#######################################################################
-!#######################################################################
-!#######################################################################
-
- FUNCTION HSL_f(t) RESULT(f)
-!=======================================================================
-! Integrand for the high-speed limit approach (Lin-Liu's formulation)
-!=======================================================================
- IMPLICIT NONE
- REAL(8), INTENT(in) :: t
- REAL(8) :: f,g
-   g = sqrt(1+t*t*q2)
-   f = t**(3+r2)/g**3 * (gp1/(g+1))**r2
- END FUNCTION HSL_f
-
-!#######################################################################
-
- END MODULE green_func_p
-
-!#######################################################################

src/itm_constants.f90

@@ -1,32 +0,0 @@
-!> Module implementing the ITM physics constants
-!>
-!> Source:
-!>  based on SOLPS b2mod_constants.F
-!>    '09/12/07 xpb : source CODATA 2006 (http://www.nist.gov/)'
-!>  pulled from ets r100
-!>
-!> \author David Coster
-!>
-!> \version "$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $"
-
-module itm_constants
-
-  use itm_types
-
-  real (kind = R8), parameter :: itm_pi = 3.141592653589793238462643383280_R8
-  real (kind = R8), parameter :: itm_c = 2.99792458e8_R8         ! speed of light, m/s
-  real (kind = R8), parameter :: itm_me = 9.10938215e-31_R8      ! electron mass, kg
-  real (kind = R8), parameter :: itm_mp = 1.672621637e-27_R8     ! proton mass, kg
-  real (kind = R8), parameter :: itm_md = 3.34358320e-27_R8      ! deuteron mass, kg
-  real (kind = R8), parameter :: itm_mt = 5.00735588e-27_R8      ! triton mass, kg
-  real (kind = R8), parameter :: itm_ma = 6.64465620e-27_R8      ! alpha mass, kg
-  real (kind = R8), parameter :: itm_amu = 1.660538782e-27_R8    ! amu, kg
-  real (kind = R8), parameter :: itm_ev = 1.602176487e-19_R8
-  real (kind = R8), parameter :: itm_qe = itm_ev
-  real (kind = R8), parameter :: itm_mu0 = 4.0e-7_R8 * itm_pi
-  real (kind = R8), parameter :: itm_eps0 = 1.0_R8 / (itm_mu0 * itm_c * itm_c)
-  real (kind = R8), parameter :: itm_avogr = 6.02214179e23_R8
-  real (kind = R8), parameter :: itm_KBolt = 1.3806504e-23_R8
-  character (len=64), parameter :: itm_constants_version = '$Id: itm_constants.f90 37 2009-08-17 17:15:00Z coster $'
-
-end module itm_constants

src/itm_types.f90

@@ -1,50 +0,0 @@
-!> Module implementing the ITM basic types
-!>
-!> Source:
-!>  based on SOLPS b2mod_types.F
-!>  pulled from ets r100 and extended with input from C. Konz, T. Ribeiro & B. Scott
-!>
-!> \author David Coster
-!>
-!> \version "$Id: itm_types.f90 144 2010-10-07 09:26:24Z konz $"
-
-module itm_types
-
-  INTEGER,  PARAMETER :: ITM_I1 = SELECTED_INT_KIND (2)        ! Integer*1
-  INTEGER,  PARAMETER :: ITM_I2 = SELECTED_INT_KIND (4)        ! Integer*2
-  INTEGER,  PARAMETER :: ITM_I4 = SELECTED_INT_KIND (9)        ! Integer*4
-  INTEGER,  PARAMETER :: ITM_I8 = SELECTED_INT_KIND (18)       ! Integer*8
-  INTEGER,  PARAMETER :: R4 = SELECTED_REAL_KIND (6, 37)   ! Real*4
-  INTEGER,  PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300) ! Real*8
-
-  INTEGER,  PARAMETER :: itm_int_invalid = -999999999
-  REAL(R8), PARAMETER :: itm_r8_invalid = -9.0D40
-
-  interface itm_is_valid
-     module procedure itm_is_valid_int4, itm_is_valid_int8, itm_is_valid_real8
-  end interface
-
-contains
-
-  logical function itm_is_valid_int4(in_int)
-    implicit none
-    integer(ITM_I4) in_int
-    itm_is_valid_int4 = in_int .ne. itm_int_invalid
-    return
-  end function itm_is_valid_int4
-
-  logical function itm_is_valid_int8(in_int)
-    implicit none
-    integer(ITM_I8) in_int
-    itm_is_valid_int8 = in_int .ne. itm_int_invalid
-    return
-  end function itm_is_valid_int8
-
-  logical function itm_is_valid_real8(in_real)
-    implicit none
-    real(R8) in_real
-    itm_is_valid_real8 = abs(in_real - itm_r8_invalid) .gt. abs(itm_r8_invalid) * 1.0e-15_R8
-    return
-  end function itm_is_valid_real8
-
-end module itm_types

src/magsurf_data.f90

@@ -0,0 +1,106 @@
+module magsurf_data
+  use const_and_precisions, only : wp_
+  implicit none
+  
+  INTEGER, SAVE :: npsi, npoints !# sup mag, # punti per sup
+  INTEGER, SAVE :: njpt, nlmt
+
+  REAL(wp_), SAVE :: rarea
+
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: psicon,pstab,rhot_eq, &
+   rhotqv,bav,varea,vcurrp,vajphiav,qqv,ffc,vratja,vratjb
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: rpstab
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: vvol,rri,rbav,bmxpsi,bmnpsi
+  REAL(wp_), DIMENSION(:), ALLOCATABLE, SAVE :: tjp,tlm,ch,ch01
+
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: rcon,zcon
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cdadrhot,cdvdrhot
+  
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cvol,crri,crbav,cbmx,cbmn,carea,cfc
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: crhotq
+  REAL(wp_), DIMENSION(:,:), ALLOCATABLE, SAVE :: cratja,cratjb,cratjpl
+
+
+contains
+
+  subroutine alloc_surf_anal(ierr)
+    implicit none
+    integer, intent(out) :: ierr
+    
+    if(npsi.le.0.or.npoints.le.0) then
+      ierr = -1
+      return
+    end if
+    
+    call dealloc_surf_anal
+    allocate(psicon(npsi),rcon(npsi,npoints), &
+     zcon(npsi,npoints),stat=ierr)
+    if (ierr/=0) call dealloc_surf_anal
+  end subroutine alloc_surf_anal
+
+  subroutine dealloc_surf_anal
+    implicit none
+    if(allocated(psicon))  deallocate(psicon)
+    if(allocated(rcon))    deallocate(rcon)
+    if(allocated(zcon))    deallocate(zcon)
+  end subroutine dealloc_surf_anal
+
+
+  subroutine alloc_surfvec(ierr)
+    implicit none
+    integer, intent(out) :: ierr
+
+    if(npsi.le.0.or.npoints.le.0) then
+      ierr = -1
+      return
+    end if
+    
+    call dealloc_surfvec
+    allocate(psicon(npsi),rcon(npsi,npoints),zcon(npsi,npoints),pstab(npsi), &
+     rhot_eq(npsi),rhotqv(npsi),bav(npsi),bmxpsi(npsi),bmnpsi(npsi),varea(npsi), &
+     vvol(npsi),vcurrp(npsi),vajphiav(npsi),qqv(npsi),ffc(npsi),vratja(npsi), &
+     vratjb(npsi),rpstab(npsi),rri(npsi),rbav(npsi),cdadrhot(npsi,4), &
+     cdvdrhot(npsi,4),cbmx(npsi,4),cbmn(npsi,4),crbav(npsi,4),cvol(npsi,4), &
+     crri(npsi,4),carea(npsi,4),cfc(npsi,4),crhotq(npsi,4),cratjpl(npsi,4), &
+     cratja(npsi,4),cratjb(npsi,4),stat=ierr)
+    if (ierr/=0) call dealloc_surf_anal
+  end subroutine alloc_surfvec
+
+  subroutine dealloc_surfvec
+    implicit none
+    if(allocated(psicon))  deallocate(psicon)
+    if(allocated(rcon))  deallocate(rcon)
+    if(allocated(zcon))  deallocate(zcon)
+    if(allocated(pstab))  deallocate(pstab)
+    if(allocated(rhot_eq))  deallocate(rhot_eq)
+    if(allocated(rhotqv))  deallocate(rhotqv)
+    if(allocated(bav))  deallocate(bav)
+    if(allocated(bmxpsi))  deallocate(bmxpsi)
+    if(allocated(bmnpsi))  deallocate(bmnpsi)
+    if(allocated(varea))  deallocate(varea)
+    if(allocated(vvol))  deallocate(vvol)
+    if(allocated(vcurrp))  deallocate(vcurrp)
+    if(allocated(vajphiav))  deallocate(vajphiav)
+    if(allocated(qqv))  deallocate(qqv)
+    if(allocated(ffc))  deallocate(ffc)
+    if(allocated(vratja))  deallocate(vratja)
+    if(allocated(vratjb))  deallocate(vratjb)
+    if(allocated(rpstab))  deallocate(rpstab)
+    if(allocated(rri))  deallocate(rri)
+    if(allocated(rbav))  deallocate(rbav)
+    if(allocated(cdadrhot))  deallocate(cdadrhot)
+    if(allocated(cdvdrhot))  deallocate(cdvdrhot)
+    if(allocated(cbmx))  deallocate(cbmx)
+    if(allocated(cbmn))  deallocate(cbmn)
+    if(allocated(crbav))  deallocate(crbav)
+    if(allocated(cvol))  deallocate(cvol)
+    if(allocated(crri))  deallocate(crri)
+    if(allocated(carea))  deallocate(carea)
+    if(allocated(cfc))  deallocate(cfc)
+    if(allocated(crhotq))  deallocate(crhotq)
+    if(allocated(cratjpl))  deallocate(cratjpl)
+    if(allocated(cratja))  deallocate(cratja)
+    if(allocated(cratjb))  deallocate(cratjb)
+  end subroutine dealloc_surfvec
+
+end module magsurf_data

src/math.f90

@@ -0,0 +1,125 @@
+module math
+
+  use const_and_precisions, only : wp_, zero, one
+  implicit none
+
+contains
+
+  function catand(z)
+!***begin prologue  catan
+!***purpose  compute the complex arc tangent.
+!***library   slatec (fnlib)
+!***category  c4a
+!***type      complex (catan-c)
+!***keywords  arc tangent, elementary functions, fnlib, trigonometric
+!***author  fullerton, w., (lanl)
+!***description
+!
+! catan(z) calculates the complex trigonometric arc tangent of z.
+! the result is in units of radians, and the real part is in the first
+! or fourth quadrant.
+!
+!***references  (none)
+!***routines called  (none)
+!***revision history  (yymmdd)
+!   770801  date written
+!   890531  changed all specific intrinsics to generic.  (wrb)
+!   890531  revision date from version 3.2
+!   891214  prologue converted to version 4.0 format.  (bab)
+!   900315  calls to xerror changed to calls to xermsg.  (thj)
+!   900326  removed duplicate information from description section.
+!       (wrb)
+!***end prologue  catan
+    use const_and_precisions, only : comp_eps, pi2=>pihalf, czero, cunit
+    implicit none
+    complex(wp_) :: catand
+    complex(wp_), intent(in) :: z
+    complex(wp_) :: z2
+    real(wp_) :: r,x,y,r2,xans,yans,twoi
+    integer :: i
+    logical, save :: first=.true.
+    integer, save :: nterms
+    real(wp_), save :: rmin, rmax, sqeps
+!***first executable statement  catan
+    if (first) then
+! nterms = log(eps)/log(rbnd) where rbnd = 0.1
+      nterms = int(-0.4343_wp_*log(0.5_wp_*comp_eps) + 1.0_wp_)
+      sqeps = sqrt(comp_eps)
+      rmin = sqrt (1.5_wp_*comp_eps)
+      rmax = 2.0_wp_/comp_eps
+    endif
+    first = .false.
+!
+    r = abs(z)
+    if (r<=0.1_wp_) then
+!
+      catand = z
+      if (r<rmin) return
+!
+      catand = czero
+      z2 = z*z
+      do i=1,nterms
+        twoi = 2*(nterms-i) + 1
+        catand = 1.0_wp_/twoi - z2*catand
+      end do
+      catand = z*catand
+!
+    else if (r<=rmax) then
+      x = real(z)
+      y = aimag(z)
+      r2 = r*r
+      if (r2==one.and.x==zero) print*,'catand, z is +i or -i'
+      if (abs(r2-one)<=sqeps) then
+        if (abs(cunit+z*z) < sqeps) &
+         print*,'catand, answer lt half precision, z**2 close to -1'
+!
+      end if
+      xans = 0.5_wp_*atan2(2.0_wp_*x, one)
+      yans = 0.25_wp_*log((r2+2.0_wp_*y+one)/(r2-2.0_wp_*y+one))
+      catand = cmplx(xans, yans, wp_)
+!
+    else
+      catand = cmplx(pi2, zero, wp_)
+      if (real(z)<zero) catand = cmplx(-pi2, zero, wp_)
+    end if
+  end function catand
+
+  function fact(k)
+    implicit none
+    integer, intent(in) :: k
+    real(wp_) :: fact
+    integer :: i
+!   Factorial function
+    fact=zero
+    if(k<0) return
+    fact=one
+    if(k==0) return
+    do i=1,k
+      fact=fact*i
+    end do
+  end function fact
+
+  function gamm(xx)
+    implicit none
+    real(wp_) :: gamm
+    real(wp_), intent(in) :: xx
+!   Returns the value Gamma(xx) for xx > 0.
+    INTEGER :: j
+    real(wp_) :: ser,tmp,x,y
+    real(wp_), parameter :: stp=2.5066282746310005_wp_
+    real(wp_), dimension(6), parameter :: cof=(/76.18009172947146_wp_, &
+    -86.50532032941677_wp_,24.01409824083091_wp_,-1.231739572450155_wp_, &
+    .1208650973866179e-2_wp_,-.5395239384953e-5_wp_/)
+    x=xx
+    y=x
+    tmp=x+5.5_wp_
+    tmp=(x+0.5_wp_)*log(tmp)-tmp
+    ser=1.000000000190015_wp_
+    do  j=1,6
+      y=y+1._wp_
+      ser=ser+cof(j)/y
+    end do
+    gamm=exp(tmp)*(stp*ser/x)
+  end function gamm
+
+end module math

src/numint.f90

@@ -0,0 +1,257 @@
+module numint
+
+  use const_and_precisions, only : wp_, zero, one
+  implicit none
+
+contains
+
+  subroutine simpson (n,h,fi,s)     
+!   subroutine for integration over f(x) with the simpson rule.  fi:
+!   integrand f(x); h: interval; s: integral.  copyright (c) tao pang 1997.
+    implicit none
+    integer, intent(in) ::  n
+    real(wp_), intent(in) :: h
+    real(wp_), dimension(n), intent(in) :: fi
+    real(wp_), intent(out) ::  s
+    integer  :: i
+    real(wp_) :: s0,s1,s2
+
+    s  = zero
+    s0 = zero
+    s1 = zero
+    s2 = zero
+    do i = 2, n-1, 2
+      s1 = s1+fi(i-1)
+      s0 = s0+fi(i)
+      s2 = s2+fi(i+1)
+    end do
+    s = h*(s1+4.0_wp_*s0+s2)/3.0_wp_
+! if n is even, add the last slice separately     
+    if (mod(n,2).eq.0) s = s+h*(5.0_wp_*fi(n)+8.0_wp_*fi(n-1)-fi(n-2))/12.0_wp_
+  end subroutine simpson
+
+  subroutine trapezoid(n,xi,fi,s)
+!   subroutine for integration with the trapezoidal rule.
+!   fi: integrand f(x); xi: abscissa x;
+!   s: integral Int_{xi(1)}^{xi(n)} f(x)dx
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), dimension(n), intent(in) :: xi,fi
+    real(wp_), intent(out) :: s
+    integer :: i
+
+    s  = zero
+    do i = 1, n-1
+      s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
+    end do
+    s = 0.5_wp_*s
+  end subroutine trapezoid
+
+  subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
+    implicit none
+    real(wp_), intent(in) :: a, b, abserr, relerr
+    real(wp_), intent(out) :: result, errest, flag
+    integer, intent(out) :: nofun
+!
+!   estimate the integral of fun(x) from a to b
+!   to a user provided tolerance.
+!   an automatic adaptive routine based on
+!   the 8-panel newton-cotes rule.
+!
+!   input ..
+!
+!   fun     the name of the integrand function subprogram fun(x).
+!   a       the lower limit of integration.
+!   b       the upper limit of integration.(b may be less than a.)
+!   relerr  a relative error tolerance. (should be non-negative)
+!   abserr  an absolute error tolerance. (should be non-negative)
+!
+!   output ..
+!
+!   result  an approximation to the integral hopefully satisfying the
+!           least stringent of the two error tolerances.
+!   errest  an estimate of the magnitude of the actual error.
+!   nofun   the number of function values used in calculation of result.
+!   flag    a reliability indicator.  if flag is zero, then result
+!           probably satisfies the error tolerance.  if flag is
+!           xxx.yyy , then  xxx = the number of intervals which have
+!           not converged and  0.yyy = the fraction of the interval
+!           left to do when the limit on  nofun  was approached.
+!
+    real(wp_) :: w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
+    real(wp_) :: qprev,qnow,qdiff,qleft,esterr,tolerr
+    real(wp_), dimension(31) :: qright
+    real(wp_), dimension(16) :: f,x
+    real(wp_), dimension(8,30) :: fsave,xsave
+    integer :: levmin,levmax,levout,nomax,nofin,lev,nim,i,j
+
+    interface
+       function fun(x)
+         use const_and_precisions, only : wp_
+         implicit none
+         real(wp_), intent(in) :: x
+         real(wp_) :: fun
+       end function fun
+    end interface
+!
+!   ***   stage 1 ***   general initialization
+!   set constants.
+!
+    levmin = 1
+    levmax = 30
+    levout = 6
+    nomax = 5000
+    nofin = nomax - 8*(levmax-levout+2**(levout+1))
+!
+!   trouble when nofun reaches nofin
+!
+    w0 =   3956.0_wp_ / 14175.0_wp_
+    w1 =  23552.0_wp_ / 14175.0_wp_
+    w2 =  -3712.0_wp_ / 14175.0_wp_
+    w3 =  41984.0_wp_ / 14175.0_wp_
+    w4 = -18160.0_wp_ / 14175.0_wp_
+!
+!   initialize running sums to zero.
+!
+    flag = zero
+    result = zero
+    cor11  = zero
+    errest = zero
+    area   = zero
+    nofun = 0
+    if (a .eq. b) return
+!
+!   ***   stage 2 ***   initialization for first interval
+!
+    lev = 0
+    nim = 1
+    x0 = a
+    x(16) = b
+    qprev  = zero
+    f0 = fun(x0)
+    stone = (b - a) / 16.0_wp_
+    x(8)  =  (x0  + x(16)) / 2.0_wp_
+    x(4)  =  (x0  + x(8))  / 2.0_wp_
+    x(12) =  (x(8)  + x(16)) / 2.0_wp_
+    x(2)  =  (x0  + x(4))  / 2.0_wp_
+    x(6)  =  (x(4)  + x(8))  / 2.0_wp_
+    x(10) =  (x(8)  + x(12)) / 2.0_wp_
+    x(14) =  (x(12) + x(16)) / 2.0_wp_
+    do j = 2, 16, 2
+       f(j) = fun(x(j))
+    end do
+    nofun = 9
+!
+!   ***   stage 3 ***   central calculation
+!   requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16.
+!   calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area.
+!
+    do
+      do
+        x(1) = (x0 + x(2)) / 2.0_wp_
+        f(1) = fun(x(1))
+        do j = 3, 15, 2
+          x(j) = (x(j-1) + x(j+1)) / 2.0_wp_
+          f(j) = fun(x(j))
+        end do
+        nofun = nofun + 8
+        step = (x(16) - x0) / 16.0_wp_
+        qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) &
+           + w3*(f(3)+f(5)) + w4*f(4)) * step
+        qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) &
+           + w3*(f(11)+f(13)) + w4*f(12)) * step
+        qnow = qleft + qright(lev+1)
+        qdiff = qnow - qprev
+        area = area + qdiff
+!
+!   ***   stage 4 *** interval convergence test
+!
+        esterr = abs(qdiff) / 1023.0_wp_
+        tolerr = max(abserr,relerr*abs(area)) * (step/stone)
+        if (lev .ge. levmin) then
+!
+!   ***   stage 6   ***   trouble section
+!   number of function values is about to exceed limit.
+!
+          if (lev .ge. levmax) then
+!
+!   current level is levmax.
+!
+            flag = flag + one
+            exit
+          end if
+          if (nofun .gt. nofin) then
+            nofin = 2*nofin
+            levmax = levout
+            flag = flag + (b - x0) / (b - a)
+            exit
+          end if
+          if (esterr .le. tolerr) exit
+        end if
+!
+!   ***   stage 5   ***   no convergence
+!   locate next interval.
+!
+        nim = 2*nim
+        lev = lev+1
+!
+!   store right hand elements for future use.
+!
+        do i = 1, 8
+          fsave(i,lev) = f(i+8)
+          xsave(i,lev) = x(i+8)
+        end do
+!
+!   assemble left hand elements for immediate use.
+!
+        qprev = qleft
+        do i = 1, 8
+          j = -i
+          f(2*j+18) = f(j+9)
+          x(2*j+18) = x(j+9)
+        end do
+      end do
+!
+!   ***   stage 7   ***   interval converged
+!   add contributions into running sums.
+!
+      result = result + qnow
+      errest = errest + esterr
+      cor11  = cor11  + qdiff / 1023.0_wp_
+!
+!   locate next interval.
+!
+      do
+        if (nim .eq. 2*(nim/2)) exit
+        nim = nim/2
+        lev = lev-1
+      end do
+      nim = nim + 1
+      if (lev .le. 0) exit
+!
+!   assemble elements required for the next interval.
+!
+      qprev = qright(lev)
+      x0 = x(16)
+      f0 = f(16)
+      do i = 1, 8
+        f(2*i) = fsave(i,lev)
+        x(2*i) = xsave(i,lev)
+      end do
+    end do
+!
+!   ***   stage 8   ***   finalize and return
+!
+    result = result + cor11
+!
+!   make sure errest not less than roundoff level.
+!
+    if (errest .eq. zero) return
+    do
+      temp = abs(result) + errest
+      if (temp .ne. abs(result)) return
+      errest = 2.0_wp_*errest
+    end do
+  end subroutine quanc8
+
+end module numint

src/reflections.f90

@@ -1,21 +1,27 @@
 module reflections
+  use const_and_precisions, only : wp_, comp_tiny, comp_eps, comp_huge, zero, one
   implicit none
+
+  ! === 1D array limiter Rlim_i, Zlim_i
+  integer, public, save :: nlim
+  real(wp_), public, dimension(:), allocatable, save :: rlim,zlim
+
   private
-  integer, parameter :: r8=selected_real_kind(15,300)
-  real(r8), parameter :: tinyr8=tiny(1._r8)
   public :: reflect,inters_linewall,inside
   public :: linecone_coord,interssegm_coord,interssegm
+  public :: alloc_lim,wall_refl
+
 contains
 
 subroutine reflect(ki,nsurf,ko)
   implicit none
-  real(r8), intent(in), dimension(3) :: ki
-  real(r8), intent(in), dimension(3) :: nsurf
-  real(r8), intent(out), dimension(3) :: ko
-  real(r8) :: twokn,norm2
+  real(wp_), intent(in), dimension(3) :: ki
+  real(wp_), intent(in), dimension(3) :: nsurf
+  real(wp_), intent(out), dimension(3) :: ko
+  real(wp_) :: twokn,norm2
   norm2 = dot_product(nsurf,nsurf)
-  if (norm2>0.0_r8) then
-    twokn = 2.0_r8*dot_product(ki,nsurf)/norm2
+  if (norm2>zero) then
+    twokn = 2.0_wp_*dot_product(ki,nsurf)/norm2
     ko=ki-twokn*nsurf
   else
     ko=ki
@@ -24,30 +30,38 @@ end subroutine reflect
 
 subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   implicit none
-  real(r8), intent(in), dimension(3) :: xv,kv
+  real(wp_), intent(in), dimension(3) :: xv,kv
   integer, intent(in) :: nw
-  real(r8), dimension(nw), intent(in) :: rw,zw
-  real(r8), intent(out) :: sint
-  real(r8), dimension(3), intent(out) :: normw
-  integer :: i,j,ni,iint
-  real(r8), dimension(2) :: si,ti
-  real(r8) :: drw,dzw,xint,yint,rint,l,kxy
-  real(r8) :: tol
-  tol=sqrt(epsilon(1.0_r8))
-  sint=huge(sint)
+  real(wp_), dimension(nw), intent(in) :: rw,zw
+  real(wp_), intent(out) :: sint
+  real(wp_), dimension(3), intent(out) :: normw
+  integer :: i,j,ni,iint,nneg
+  real(wp_), dimension(2) :: si,ti
+  real(wp_) :: drw,dzw,xint,yint,rint,l,kxy
+  real(wp_) :: tol
+  tol=sqrt(comp_eps)
+  sint=comp_huge
   iint=0
-  normw=0.0_r8
+  normw=zero
   do i=1,nw-1
     !search intersections with i-th wall segment
     call linecone_coord(xv,kv,rw(i:i+1),zw(i:i+1),si,ti,ni)
-    do while (ni>0 .and. si(1)<=tol)
-      !remove solutions with s<=0
-      ni = ni-1
-      si(1) = si(2)
-      ti(1) = ti(2)
-    end do
+    !discard solutions with s<=0
+    nneg=0
     do j=1,ni
-      if ((si(j)<sint .or. iint==0) .and. ti(j)>=0._r8 .and. ti(j)<=1._r8) then
+      if (si(j)<=tol) then
+        nneg=j
+      else
+        exit
+      end if
+    end do
+!    do while (ni>0 .and. si(1)<=tol)
+!      ni = ni-1
+!      si(1) = si(2)   ???
+!      ti(1) = ti(2)   ???
+!    end do
+    do j=nneg+1,ni
+      if ((si(j)<sint .or. iint==0) .and. ti(j)>=zero .and. ti(j)<=one) then
         !check intersection is in r,z range and keep the closest
         sint = si(j)
         iint = i
@@ -64,7 +78,7 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
   l = sqrt(drw**2+dzw**2)
   kxy = sqrt(kv(1)**2+kv(2)**2)
   normw(3) = -drw/l
-  if (rint>0.0_r8) then
+  if (rint>zero) then
     normw(1) = xint/rint*dzw/l
     normw(2) = yint/rint*dzw/l
   else
@@ -72,17 +86,18 @@ subroutine inters_linewall(xv,kv,rw,zw,nw,sint,normw)
     normw(2) = kv(2)/kxy*dzw/l
   end if
   !reverse normal if k.n>0
-  if (dot_product(normw,kv)>0.0_r8) normw=-normw
+  if (dot_product(normw,kv)>zero) normw=-normw
 end subroutine inters_linewall
 
 subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
+  use utils, only : bubble
   implicit none
-  real(r8), intent(in), dimension(3) :: xv,kv
-  real(r8), intent(in), dimension(2) :: rs,zs
-  real(r8), dimension(2), intent(out) :: s,t
+  real(wp_), intent(in), dimension(3) :: xv,kv
+  real(wp_), intent(in), dimension(2) :: rs,zs
+  real(wp_), dimension(2), intent(out) :: s,t
   integer, intent(out) :: n
-  real(r8) :: x0,y0,z0,kx,ky,kz
-  real(r8) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin
+  real(wp_) :: x0,y0,z0,kx,ky,kz
+  real(wp_) :: dr,dz,r,a,bhalf,c,delta,tvertex,zvertex,srmin,rmin,zrmin
   x0=xv(1)
   y0=xv(2)
   z0=xv(3)
@@ -93,9 +108,9 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
   dz = zs(2)-zs(1)
   s = 0
   t = 0
-  if (abs(dz)<tinyr8) then
+  if (abs(dz)<comp_tiny) then
     !surface in horizontal plane
-    if (abs(kz)<tinyr8 .or. abs(dr)<tinyr8) then
+    if (abs(kz)<comp_tiny .or. abs(dr)<comp_tiny) then
       n = 0
     else
       s(1) = (zs(1)-z0)/kz
@@ -107,9 +122,9 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
     a = (kx**2+ky**2) - (dr/dz*kz)**2
     bhalf = -dr/dz*kz*rs(1) + (kx*x0 + ky*y0) - (dr/dz)**2*kz*(z0-zs(1))
     c = (x0**2+y0**2) - (rs(1) + dr/dz*(z0-zs(1)))**2
-    if (abs(a)<tinyr8) then
+    if (abs(a)<comp_tiny) then
       !line parallel to cone generator
-      if (abs(dr)<tinyr8) then
+      if (abs(dr)<comp_tiny) then
         !cylinder and vertical line
         n = 0
       else
@@ -118,14 +133,14 @@ subroutine linecone_coord(xv,kv,rs,zs,s,t,n)
         srmin = -(kx*x0 + ky*y0)/(kx**2+ky**2)
         rmin = sqrt((x0+srmin*kx)**2+(y0+srmin*ky)**2)
         zrmin = z0 + srmin*kz
-        if (rmin<tinyr8 .and. abs(zrmin-zvertex)<tinyr8) then
+        if (rmin<comp_tiny .and. abs(zrmin-zvertex)<comp_tiny) then
           !line passing by cone vertex
           !s(1) = srmin
           !t(1) = tvertex
           !n = 1
           n = 0
         else
-          s(1) = -0.5_r8*c/bhalf
+          s(1) = -0.5_wp_*c/bhalf
           t(1) = (kz*s(1)+(z0-zs(1)))/dz
           n = 1
         end if
@@ -147,18 +162,18 @@ end subroutine linecone_coord
 
 subroutine interssegm_coord(xa,ya,xb,yb,s,t,ierr)
   implicit none
-  real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
-  real(r8), intent(out) :: s,t
+  real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
+  real(wp_), intent(out) :: s,t
   integer, intent(out) :: ierr
-  real(r8) :: crossprod,dxa,dya,dxb,dyb
+  real(wp_) :: crossprod,dxa,dya,dxb,dyb
   dxa = xa(2)-xa(1)
   dya = ya(2)-ya(1)
   dxb = xb(2)-xb(1)
   dyb = yb(2)-yb(1)
   crossprod = dxb*dya - dxa*dyb
-  if (abs(crossprod)<tiny(crossprod)) then
-    s = 0.0_r8
-    t = 0.0_r8
+  if (abs(crossprod)<comp_tiny) then
+    s = zero
+    t = zero
     ierr = 1
   else
     s = (dyb*(xa(1)-xb(1)) - dxb*(ya(1)-yb(1)))/crossprod
@@ -169,25 +184,26 @@ end subroutine interssegm_coord
 
 function interssegm(xa,ya,xb,yb)
   implicit none
-  real(r8), dimension(2), intent(in) :: xa,ya,xb,yb
+  real(wp_), dimension(2), intent(in) :: xa,ya,xb,yb
   logical :: interssegm
-  real(r8) :: s,t
+  real(wp_) :: s,t
   integer :: ierr
   interssegm = .false.
   call interssegm_coord(xa,ya,xb,yb,s,t,ierr)
-  if (ierr==0 .and. s>=0._r8 .and. s<=1._r8 .and. &
-                    t>=0._r8 .and. t<=1._r8) interssegm = .true.
+  if (ierr==0 .and. s>=zero .and. s<=one .and. &
+                    t>=zero .and. t<=one) interssegm = .true.
 end function interssegm
 
 function inside(xc,yc,n,x,y)
+  use utils, only : locatef, locate_unord, intlinf, bubble
   implicit none
   integer, intent(in) :: n
-  real(r8), dimension(n), intent(in) :: xc,yc
-  real(r8), intent(in) :: x,y
+  real(wp_), dimension(n), intent(in) :: xc,yc
+  real(wp_), intent(in) :: x,y
   logical :: inside
   integer, dimension(n) :: jint
-  real(r8), dimension(n) :: xint
-  real(r8), dimension(n+1) :: xclosed,yclosed
+  real(wp_), dimension(n) :: xint
+  real(wp_), dimension(n+1) :: xclosed,yclosed
   integer :: i,nj
   xclosed(1:n)=xc(1:n)
   yclosed(1:n)=yc(1:n)
@@ -197,92 +213,107 @@ function inside(xc,yc,n,x,y)
   inside=.false.
   if (nj==0) return
   do i=1,nj
-    xint(i)=intlin(yclosed(jint(i)),xclosed(jint(i)), &
+    xint(i)=intlinf(yclosed(jint(i)),xclosed(jint(i)), &
                    yclosed(jint(i)+1),xclosed(jint(i)+1),y)
   end do
   call bubble(xint,nj)
-  inside=(mod(locate(xint,nj,x),2)==1)
+  inside=(mod(locatef(xint,nj,x),2)==1)
 end function inside
 
-function intlin(x1,y1,x2,y2,x) result(y)
-  !linear interpolation
-  !must be x1 != x2
+subroutine alloc_lim(ier)
   implicit none
-  real(r8),intent(in) :: x1,y1,x2,y2,x
-  real(r8) :: y
-  real(r8) :: a
-  a=(x2-x)/(x2-x1)
-  y=a*y1+(1._r8-a)*y2
-end function intlin
-
-subroutine locate_unord(a,n,x,j,m,nj)
-  implicit none
-  integer, intent(in) :: n,m
-  integer, intent(out) :: nj
-  real(r8), dimension(n), intent(in) :: a
-  real(r8), intent(in) :: x
-  integer, dimension(m), intent(inout) :: j
-  integer :: i
-  nj=0
-  do i=1,n-1
-    if (x>a(i).neqv.x>a(i+1)) then
-      nj=nj+1
-      if (nj<=m) j(nj)=i
-    end if
-  end do
-end subroutine locate_unord
-
-function locate(a,n,x) result(j)
-  !Given an array a(n), and a value x, with a(n) monotonic, either
-  !increasing or decreasing, returns a value j such that
-  !a(j) < x <= a(j+1) for a increasing, and such that
-  !a(j+1) < x <= a(j) for a decreasing.
-  !j=0 or j=n  indicate that x is out of range  (Numerical Recipes)
-  implicit none
-  integer, intent(in) :: n
-  real(r8), dimension(n), intent(in) :: a
-  real(r8), intent(in) :: x
-  integer :: j
-  integer :: jl,ju,jm
-  logical :: incr
-  jl=0
-  ju=n+1
-  incr=a(n)>a(1)
-  do while ((ju-jl)>1)
-    jm=(ju+jl)/2
-    if(incr.eqv.(x>a(jm))) then
-      jl=jm
-    else
-      ju=jm
-    endif
-  end do
-  j=jl
-end function locate
-
-subroutine order(p,q)
-  !returns p,q in ascending order
-  implicit none
-  real(r8), intent(inout) :: p,q
-  real(r8) :: temp
-  if (p>q) then
-    temp=p
-    p=q
-    q=temp
+  integer, intent(out) :: ier
+  
+  if(nlim.lt.0) then
+    ier = -1
+    return
   end if
-end subroutine order
+  
+  call dealloc_lim
+  allocate(rlim(nlim),zlim(nlim), &
+   stat=ier)
+  if (ier/=0) call dealloc_lim
+end subroutine alloc_lim
 
-subroutine bubble(a,n)
-  !bubble sorting of array a
+subroutine dealloc_lim
   implicit none
-  integer, intent(in) :: n
-  real(r8), dimension(n), intent(inout) :: a
-  integer :: i, j
-  do i=1,n
-    do j=n,i+1,-1
-      call order(a(j-1), a(j))
-    end do
-  end do
-end subroutine bubble
+  if(allocated(rlim))   deallocate(rlim)
+  if(allocated(zlim))   deallocate(zlim)
+end subroutine dealloc_lim
+
+subroutine wall_refl(xv,anv,ext,eyt,xvrfl,anvrfl,extr,eytr,walln,irfl)
+  implicit none
+! arguments
+  integer :: irfl
+  real(wp_), dimension(3) :: xv,anv,xvrfl,anvrfl,walln
+  complex(wp_) :: ext,eyt,extr,eytr
+! local variables
+  real(wp_) :: smax,rrm,zzm
+  real(wp_), dimension(3) :: anv0,vv1,vv2,vv3
+  complex(wp_) :: eztr
+  complex(wp_), dimension(3) :: evin,evrfl 
+!
+  anv0=anv/sqrt(anv(1)**2+anv(2)**2+anv(3)**2)
+  rrm=1.0e-2_wp_*sqrt(xv(1)**2+xv(2)**2)
+  zzm=1.0e-2_wp_*xv(3)
+!
+! computation of reflection coordinates and normal to the wall
+  call inters_linewall(xv/1.0e2_wp_,anv0,rlim(1:nlim),zlim(1:nlim), &
+                       nlim,smax,walln)
+  smax=smax*1.0e2_wp_
+  xvrfl=xv+smax*anv0
+  irfl=1
+  if (.not.inside(rlim,zlim,nlim,rrm,zzm)) then
+    ! first wall interface is outside-inside
+    if (dot_product(walln,walln)<tiny(walln)) then
+      ! wall never hit
+      xvrfl=xv
+      anvrfl=anv0
+      extr=ext
+      eytr=eyt
+      irfl=0
+      return
+    end if
+    ! search second wall interface (inside-outside)
+    call inters_linewall(xvrfl/1.0e2_wp_,anv0,rlim(1:nlim), &
+                         zlim(1:nlim),nlim,smax,walln)
+    smax=smax*1.0e2_wp_
+    xvrfl=xvrfl+smax*anv0
+    irfl=2
+  end if
+!
+! rotation matrix from local to lab frame
+  vv1(1)=anv0(2)
+  vv1(2)=-anv0(1)
+  vv1(3)=0.0_wp_
+  vv2(1)=anv0(1)*anv0(3)
+  vv2(2)=anv0(2)*anv0(3)
+  vv2(3)=-anv0(1)*anv0(1)-anv0(2)*anv0(2)
+  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
+  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
+  vv3=anv0
+!
+  evin=ext*vv1+eyt*vv2
+! wave vector and electric field after reflection in lab frame
+  anvrfl=anv0-2.0_wp_* &
+         (anv0(1)*walln(1)+anv0(2)*walln(2)+anv0(3)*walln(3))*walln
+  evrfl=-evin+2.0_wp_* &
+         (evin(1)*walln(1)+evin(2)*walln(2)+evin(3)*walln(3))*walln
+!
+  vv1(1)=anvrfl(2)
+  vv1(2)=-anvrfl(1)
+  vv1(3)=0.0_wp_
+  vv2(1)=anvrfl(1)*anvrfl(3)
+  vv2(2)=anvrfl(2)*anvrfl(3)
+  vv2(3)=-anvrfl(1)*anvrfl(1)-anvrfl(2)*anvrfl(2)
+  vv1=vv1/sqrt(vv1(1)**2+vv1(2)**2+vv1(3)**2)
+  vv2=vv2/sqrt(vv2(1)**2+vv2(2)**2+vv2(3)**2)
+  vv3=anvrfl/sqrt(anvrfl(1)**2+anvrfl(2)**2+anvrfl(3)**2)
+!
+  extr=dot_product(vv1,evrfl)
+  eytr=dot_product(vv2,evrfl)
+  eztr=dot_product(vv3,evrfl)
+end
 
 end module reflections
 

src/simplespline.f90

@@ -0,0 +1,273 @@
+module simplespline
+
+  use const_and_precisions, only : wp_
+  implicit none
+
+contains
+
+function spli(cspli,n,k,dx)
+  implicit none
+  integer, intent(in) :: n, k
+  real(wp_), intent(in) :: cspli(n,4), dx
+  real(wp_) :: spli
+  spli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4)))
+end function spli
+
+function splid(cspli,n,k,dx)
+  implicit none
+  integer, intent(in) :: n, k
+  real(wp_), intent(in) :: cspli(n,4), dx
+  real(wp_) :: splid
+  splid=cspli(k,2)+dx*(2.0_wp_*cspli(k,3)+3.0_wp_*dx*cspli(k,4))
+end function splid
+
+subroutine difcs(x,y,n,iopt,c,ier)
+  implicit none
+  integer, intent(in) :: n, iopt
+  real(wp_), intent(in) :: x(n), y(n)
+  real(wp_), intent(inout) :: c(n*4)
+  integer :: ier
+  integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
+  real(wp_) :: xb,xc,ya,yb,h,a,r,dya,dyb,dy2
+  jmp  =1
+  if (n <= 1) return
+!
+! initialization
+!
+  xc   =x(1)
+  yb   =y(1)
+  h    =0.0_wp_
+  a    =0.0_wp_
+  r    =0.0_wp_
+  dyb  =0.0_wp_
+!
+! iol=0 - given derivative at first point
+! ioh=0 - given derivative at last point
+!
+  iol  =iopt-1
+  ioh  =iopt-2
+  if (ioh == 1) then
+    iol  =0
+    ioh  =0
+  end if
+  dy2  =c(2)
+!
+! form the system of linear equations
+! and eliminate subsequentially
+!
+  j    =1
+  do i=1,n
+    j2   =n+i
+    j3   =j2+n
+    a    =h*(2.0_wp_-a)
+    dya  =dyb+h*r
+    if (i>=n) then
+!
+! set derivative dy2 at last point
+!
+      dyb  =dy2
+      h    =0.0_wp_
+      if (ioh/=0) then
+        dyb  =dya
+        goto 13
+      end if
+    else
+      j    =j+jmp
+      xb   =xc
+      xc   =x(j)
+      h    =xc-xb
+!
+! ii=0 - increasing abscissae
+! ii=1 - decreasing abscissae
+!
+      ii   =0
+      if (h==0) return
+      if (h<0) ii   =1
+      ya   =yb
+      yb   =y(j)
+      dyb  =(yb-ya)/h
+      if (i<=1) then
+        j1   =ii
+        if (iol/=0) goto 13
+        dya  =c(1)
+      end if
+    end if
+    if (j1-ii /= 0) return
+    a    =1.0_wp_/(h+h+a)
+    13 continue
+    r    =a*(dyb-dya)
+    c(j3)=r
+    a    =h*a
+    c(j2)=a
+    c(i) =dyb
+  end do
+!
+! back substitution of the system of linear equations
+! and computation of the other coefficients
+!
+  a    =1.0_wp_
+  j1   =j3+n+ii-ii*n
+  i    =n
+  do iol=1,n
+    xb   =x(j)
+    h    =xc-xb
+    xc   =xb
+    a    =a+h
+    yb   =r
+    r    =c(j3)-r*c(j2)
+    ya   =r+r
+    c(j3)=ya+r
+    c(j2)=c(i)-h*(ya+yb)
+    c(j1)=(yb-r)/a
+    c(i) =y(j)
+    a    =0.0_wp_
+    j    =j-jmp
+    i    =i-1
+    j2   =j2-1
+    j3   =j3-1
+    j1   =j3+n+ii
+  end do
+  ier  =0
+end subroutine difcs
+
+subroutine difcsn(xx,yy,nmx,n,iopt,cc,ier)
+!
+! same as difcs but with dimension(xx,yy) = nmx > n
+!
+  implicit none
+  integer, intent(in) :: nmx, n, iopt
+  real(wp_), intent(in) :: xx(nmx), yy(nmx)
+  real(wp_), intent(inout) :: cc(nmx,4)
+  integer :: ier
+  integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
+  real(wp_) :: x(n),y(n),c(n*4),xb,xc,ya,yb,h,a,r,dya,dyb,dy2
+!
+  do i=1,n
+    x(i)=xx(i)
+    y(i)=yy(i)
+  end do
+  ii=0
+  do j=1,4
+    do i=1,n
+      ii=ii+1
+      c(ii)=cc(i,j)
+    end do
+  end do
+!
+  jmp  =1
+  if (n>1) then
+!
+! initialization
+!
+    xc   =x(1)
+    yb   =y(1)
+    h    =0.0_wp_
+    a    =0.0_wp_
+    r    =0.0_wp_
+    dyb  =0.0_wp_
+!
+! iol=0 - given derivative at first point
+! ioh=0 - given derivative at last point
+!
+    iol  =iopt-1
+    ioh  =iopt-2
+    if (ioh==1) then
+      iol  =0
+      ioh  =0
+    end if
+    dy2  =c(2)
+!
+! form the system of linear equations
+! and eliminate subsequentially
+!
+    j    =1
+    do i=1,n
+      j2   =n+i
+      j3   =j2+n
+      a    =h*(2.0_wp_-a)
+      dya  =dyb+h*r
+      if (i>=n) then
+!
+! set derivative dy2 at last point
+!
+        dyb  =dy2
+        h    =0.0_wp_
+        if (ioh/=0) then
+          dyb  =dya
+          goto 13
+        end if
+      else
+        j    =j+jmp
+        xb   =xc
+        xc   =x(j)
+        h    =xc-xb
+!
+! ii=0 - increasing abscissae
+! ii=1 - decreasing abscissae
+!
+        ii   =0
+        if (h==0) goto 16
+        if (h<0) ii   =1
+        ya   =yb
+        yb   =y(j)
+        dyb  =(yb-ya)/h
+        if (i<=1) then
+          j1   =ii
+          if (iol/=0) goto 13
+          dya  =c(1)
+        end if
+      end if
+      if (j1/=ii) goto 16
+      a    =1.0_wp_/(h+h+a)
+      13 continue
+      r    =a*(dyb-dya)
+      c(j3)=r
+      a    =h*a
+      c(j2)=a
+      c(i) =dyb
+    end do
+!
+! back substitution of the system of linear equations
+! and computation of the other coefficients
+!
+    a    =1.0_wp_
+    j1   =j3+n+ii-ii*n
+    i    =n
+    do iol=1,n
+       xb   =x(j)
+       h    =xc-xb
+       xc   =xb
+       a    =a+h
+       yb   =r
+       r    =c(j3)-r*c(j2)
+       ya   =r+r
+       c(j3)=ya+r
+       c(j2)=c(i)-h*(ya+yb)
+       c(j1)=(yb-r)/a
+       c(i) =y(j)
+       a    =0.0_wp_
+       j    =j-jmp
+       i    =i-1
+       j2   =j2-1
+       j3   =j3-1
+       j1   =j3+n+ii
+    end do
+    ier  =0
+  end if
+!
+  16 continue
+  ii=0
+  do j=1,4
+    do i=1,nmx
+      if(i<=n) then
+        ii=ii+1
+        cc(i,j)=c(ii)
+       else
+         cc(i,j)=0.0_wp_
+       end if
+    end do
+  end do
+!
+end subroutine difcsn
+
+end module simplespline

src/utils.f90

@@ -0,0 +1,278 @@
+module utils
+
+  use const_and_precisions, only : wp_
+  implicit none
+
+contains
+
+  function locatef(a,n,x) result(j)
+! Given an array a(n), and a value x, with a(n) monotonic, either
+! increasing or decreasing, returns a value j such that
+! a(j) < x <= a(j+1) for a increasing, and such that
+! a(j+1) < x <= a(j) for a decreasing.
+! j=0 or j=n  indicate that x is out of range  (Numerical Recipes)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), dimension(n), intent(in) :: a
+    real(wp_), intent(in) :: x
+    integer :: j
+    integer :: jl,ju,jm
+    logical :: incr
+    jl=0
+    ju=n+1
+    incr=a(n)>a(1)
+    do while ((ju-jl)>1)
+      jm=(ju+jl)/2
+      if(incr.eqv.(x>a(jm))) then
+        jl=jm
+      else
+        ju=jm
+      endif
+    end do
+    j=jl
+  end function locatef
+
+  subroutine locate(xx,n,x,j)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), intent(in) :: xx(n), x
+    integer, intent(out) :: j
+    integer :: jl,ju,jm
+    logical :: incr
+!
+! Given an array xx(n), and a value x
+! returns a value j such that xx(j) < x < xx(j+1)
+! xx(n) must be monotonic, either increasing or decreasing.
+! j=0 or j=n  indicate that x is out of range  (Numerical Recipes)
+!
+    jl=0
+    ju=n+1
+    incr=xx(n)>xx(1)
+    do while ((ju-jl)>1)
+      jm=(ju+jl)/2
+      if(incr .eqv. (x>xx(jm))) then
+        jl=jm
+      else
+        ju=jm
+      endif
+    end do
+    j=jl
+  end subroutine locate
+
+  subroutine locatex(xx,n,n1,n2,x,j)
+    implicit none
+    integer, intent(in) :: n,n1,n2
+    real(wp_), intent(in) :: xx(n), x
+    integer, intent(out) :: j
+    integer :: jl,ju,jm
+!
+! Given an array xx(n), and a value x
+! returns a value j such that xx(j) < x < xx(j+1)
+! xx(n) must be monotonic, either increasing or decreasing.
+! j=n1-1or j=n2+1 indicate that x is out of range
+! modified from subr. locate (Numerical Recipes)
+!
+    jl=n1-1
+    ju=n2+1
+    do while ((ju-jl)>1)
+      jm=(ju+jl)/2
+      if((xx(n2)>xx(n1)) .eqv. (x>xx(jm))) then
+        jl=jm
+      else
+        ju=jm
+      endif
+    end do
+    j=jl
+  end subroutine locatex
+
+  subroutine locate_unord(a,n,x,j,m,nj)
+    implicit none
+    integer, intent(in) :: n,m
+    integer, intent(out) :: nj
+    real(wp_), dimension(n), intent(in) :: a
+    real(wp_), intent(in) :: x
+    integer, dimension(m), intent(inout) :: j
+    integer :: i
+    nj=0
+    do i=1,n-1
+      if (x>a(i).neqv.x>a(i+1)) then
+        nj=nj+1
+        if (nj<=m) j(nj)=i
+      end if
+    end do
+  end subroutine locate_unord
+
+  function intlinf(x1,y1,x2,y2,x) result(y)
+    !linear interpolation
+    !must be x1 != x2
+    use const_and_precisions, only : one
+    implicit none
+    real(wp_),intent(in) :: x1,y1,x2,y2,x
+    real(wp_) :: y
+    real(wp_) :: a
+    a=(x2-x)/(x2-x1)
+    y=a*y1+(one-a)*y2
+  end function intlinf
+
+  subroutine intlin(x1,y1,x2,y2,x,y)
+    implicit none
+    real(wp_), intent(in) :: x1,y1,x2,y2,x
+    real(wp_), intent(out) :: y
+    real(wp_) :: dx,aa,bb
+!
+! linear interpolation
+! (x1,y1) < (x,y) < (x2,y2)
+!
+    dx=x2-x1
+    aa=(x2-x)/dx
+    bb=1.0_wp_-aa
+    y=aa*y1+bb*y2
+  end subroutine intlin
+
+  subroutine vmax(x,n,xmax,imx)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), intent(in) :: x(n)
+    real(wp_), intent(out) :: xmax
+    integer, intent(out) :: imx
+    integer :: i
+
+    if (n<1) then
+      imx=0
+      return
+    end if
+    imx=1
+    xmax=x(1)
+    do i=2,n
+      if(x(i)>xmax) then
+        xmax=x(i)
+        imx=i
+      end if
+    end do
+  end subroutine vmax
+
+  subroutine vmin(x,n,xmin,imn)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), intent(in) :: x(n)
+    real(wp_), intent(out) :: xmin
+    integer, intent(out) :: imn
+    integer :: i
+
+    if (n<1) then
+      imn=0
+      return
+    end if
+    imn=1
+    xmin=x(1)
+    do i=2,n
+      if(x(i)<xmin) then
+        xmin=x(i)
+        imn=i
+      end if
+    end do
+  end subroutine vmin
+
+  subroutine vmaxmini(x,n,xmin,xmax,imn,imx)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), intent(in) :: x(n)
+    real(wp_), intent(out) :: xmin, xmax
+    integer, intent(out) :: imn, imx
+    integer :: i
+    if (n<1) then
+      imn=0
+      imx=0
+      return
+    end if
+    imn=1
+    imx=1
+    xmin=x(1)
+    xmax=x(1)
+    do i=2,n
+      if(x(i)<xmin) then
+        xmin=x(i)
+        imn=i
+      else if(x(i)>xmax) then
+        xmax=x(i)
+        imx=i
+      end if
+    end do
+  end subroutine vmaxmini
+
+  subroutine vmaxmin(x,n,xmin,xmax)
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), intent(in) :: x(n)
+    real(wp_), intent(out) :: xmin, xmax
+    integer :: i
+
+    if (n<1) then
+      return
+    end if
+    xmin=x(1)
+    xmax=x(1)
+    do i=2,n
+      if(x(i)<xmin) then
+        xmin=x(i)
+      else if(x(i)>xmax) then
+        xmax=x(i)
+      end if
+    end do
+  end subroutine vmaxmin
+
+  subroutine order(p,q)
+! returns p,q in ascending order
+    implicit none
+    real(wp_), intent(inout) :: p,q
+    real(wp_) :: temp
+    if (p>q) then
+      temp=p
+      p=q
+      q=temp
+    end if
+  end subroutine order
+
+  subroutine bubble(a,n)
+! bubble sorting of array a
+    implicit none
+    integer, intent(in) :: n
+    real(wp_), dimension(n), intent(inout) :: a
+    integer :: i, j
+    do i=1,n
+      do j=n,i+1,-1
+        call order(a(j-1), a(j))
+      end do
+    end do
+  end subroutine bubble
+
+  function get_free_unit(umin,umax) result(i)
+    implicit none
+    integer :: i
+    integer, intent(in), optional :: umin, umax
+    integer, parameter :: max_allowed = 999
+    integer :: ierr, iend
+    logical :: ex, op
+
+    if (present(umin)) then
+      i = max(0,umin)    ! start searching from unit min
+    else
+      i = 0
+    end if
+    if (present(umax)) then
+      iend = min(max(0,umax),max_allowed)
+    else
+      iend = max_allowed
+    end if
+    do
+      if (i>iend) then
+        i=-1  ! no free units found
+        exit
+      end if
+      inquire(unit=i,exist=ex,opened=op,iostat=ierr)
+      if (ierr==0.and.ex.and..not.op) exit  ! unit i exists and is not open
+      i = i + 1
+    end do
+  end function get_free_unit
+
+end module utils