spline modules added, grayl removed

Makefile

@@ -3,12 +3,11 @@ EXE=gray
 
 # Objects list
 MAINOBJ=gray.o
-OTHOBJ=conical.o const_and_precisions.o dispersion.o eierf.o \
+OTHOBJ=conical.o const_and_precisions.o dierckx.o dispersion.o eierf.o \
-    graydata_flags.o graydata_par.o graydata_anequil.o grayl.o \
+    graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
-    green_func_p.o interp_eqprof.o magsurf_data.o math.o minpack.o \
+    interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
-    numint.o quadpack.o reflections.o utils.o
+    reflections.o simplespline.o utils.o
     
-
 # Alternative search paths
 vpath %.f90 src
 vpath %.f src
@@ -26,23 +25,25 @@ $(EXE): $(MAINOBJ) $(OTHOBJ)
 	$(FC) $(FFLAGS) -o $@ $^
 
 # Dependencies on modules
-gray.o: const_and_precisions.o conical.o dispersion.o green_func_p.o \
+gray.o: const_and_precisions.o conical.o dierckx.o dispersion.o \
-        graydata_flags.o graydata_par.o graydata_anequil.o interp_eqprof.o \
+        graydata_anequil.o graydata_flags.o graydata_par.o green_func_p.o \
-        magsurf_data.o math.o minpack.o numint.o quadpack.o reflections.o \
+        interp_eqprof.o magsurf_data.o math.o minpack.o numint.o quadpack.o \
-        utils.o
+        reflections.o simplespline.o utils.o
-grayl.o: const_and_precisions.o
-green_func_p.o: const_and_precisions.o numint.o
-numint.o: const_and_precisions.o
-reflections.o: const_and_precisions.o utils.o
 conical.o: const_and_precisions.o
+dierckx.o: const_and_precisions.o
 dispersion.o: const_and_precisions.o eierf.o math.o quadpack.o
-math.o: const_and_precisions.o
+graydata_anequil.o: const_and_precisions.o
-minpack.o: const_and_precisions.o
 graydata_flags.o: const_and_precisions.o
 graydata_par.o: const_and_precisions.o
-graydata_anequil.o: const_and_precisions.o
+green_func_p.o: const_and_precisions.o numint.o
-magsurf_data.o: const_and_precisions.o
 interp_eqprof.o: const_and_precisions.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
 
 # General object compilation command

src/gray.f

@@ -1055,6 +1055,7 @@ c
       subroutine read_beams
       use graydata_flags, only : filenmbm, ibeam
       use utils, only : locate
+      use simplespline, only : difcs,spli
 
       implicit real*8(a-h,o-z)
       parameter(nstrmx=50)
@@ -1167,6 +1168,7 @@ c
       subroutine equidata
       use const_and_precisions, only : pi
       use utils, only : vmaxmini
+      use dierckx, only : curfit,splev,regrid,bispev,coeff_parder
       use graydata_flags, only : ipsinorm,sspl,ixp,icocos,neqdsk
       use graydata_par, only : sgnbphi,sgniphi,factb
       use interp_eqprof, only : nsrt,nszt,nsft,rlim,zlim,nlim,nr,nz,
@@ -1858,6 +1860,8 @@ c
       use graydata_par, only : psdbnd,factb,factt,factn
       use interp_eqprof, only : nsfd,npp,psnpp,aa,bb,cc,
      . psrad,derad,terad,zfc,tfn,cfn,ct,cz,alloc_profvec
+      use simplespline, only : difcsn
+      use dierckx, only : curfit,splev,splder
 
       implicit real*8 (a-h,o-z)
       integer, dimension(:), allocatable :: iwrkf
@@ -1968,6 +1972,7 @@ c
 c
       subroutine rhotor(nnr)
       use interp_eqprof, only : nr,psinr,qpsi,crhot,cq
+      use simplespline, only : difcsn
 
       implicit real*8(a-h,o-z)
       real*8, dimension(:), allocatable :: rhotnr
@@ -2006,6 +2011,7 @@ c
 
       function fq_eq(psi)
       use interp_eqprof, only : psinr,nr,cq
+      use simplespline, only :spli
 
       implicit real*8(a-h,o-z)
 
@@ -2019,6 +2025,7 @@ c
 
       function frhotor_eq(psi)
       use interp_eqprof, only : psinr,nr,crhot
+      use simplespline, only :spli
 
       implicit real*8(a-h,o-z)
 c
@@ -2042,6 +2049,7 @@ c      if(irt.eq.nr) irt=nr-1
       
       function frhotor_av(psi)
       use magsurf_data, only : rpstab, crhotq, npsi
+      use simplespline, only :spli
 
       implicit real*8(a-h,o-z)
 
@@ -2056,6 +2064,7 @@ c      if(ip.eq.npsi) ip=npsi-1
       end
       
       subroutine rhopol
+      use dierckx, only : curfit,splev
       implicit real*8(a-h,o-z)
       parameter(nnr=101,nrest=nnr+4)
       parameter(lwrkp=nnr*4+nrest*16)
@@ -2094,6 +2103,7 @@ c      end do
       end
       
       function frhopol(rhot)
+      use dierckx, only : splev
       implicit real*8(a-h,o-z)
       parameter(nnr=101,nrest=nnr+4)
       dimension trp(nrest),crp(nrest),rrs(1),ffspl(1)
@@ -2305,8 +2315,8 @@ c
       use magsurf_data, only : npoints
       use interp_eqprof, only : psia,psiant,psinop,nsr=>nsrt,nsz=>nszt,
      . cc=>cceq,tr,tz,rup,zup,rlw,zlw,nr
-
       use const_and_precisions, only : pi
+      use dierckx, only : profil,sproota
       implicit real*8 (a-h,o-z)
       parameter(mest=4,kspl=3)
       dimension zeroc(mest)
@@ -2370,6 +2380,8 @@ c
       use const_and_precisions, only : zero,one,pi,ccj=>mu0inv
       use magsurf_data
       use interp_eqprof, only : btrcen,btaxis,rmaxis,zmaxis,phitedge
+      use simplespline, only : difcs
+      use dierckx, only : regrid,coeff_parder
 
       implicit real*8 (a-h,o-z)
       real*8 lam
@@ -2721,6 +2733,7 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
 
       function fdadrhot(rpsi)
       use magsurf_data, only : rpstab,cdadrhot,npsi
+      use simplespline, only :spli
 
       implicit real*8(a-h,o-z)
 
@@ -2735,6 +2748,7 @@ c      if(ip.eq.npsi) ip=npsi-1
 
       function fdvdrhot(rpsi)
       use magsurf_data, only : rpstab,cdvdrhot,npsi
+      use simplespline, only :spli
 
       implicit real*8(a-h,o-z)
 
@@ -2749,6 +2763,8 @@ c      if(ip.eq.npsi) ip=npsi-1
       use graydata_anequil, only : q0,qa,alq,b0,rr0m,zr0m,rpam
       use magsurf_data
       use interp_eqprof, only : btrcen
+      use simplespline, only : difcs
+      use dierckx, only : regrid,coeff_parder
 
       implicit real*8 (a-h,o-z)
       real*8 lam
@@ -3127,6 +3143,7 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
 
 
       subroutine rhopol_an
+      use dierckx, only : curfit,splev
       use graydata_par, only : sgniphi
       use graydata_anequil, only : q0,qa,alq,b0,rr0m,rpam
       use interp_eqprof, only : psia
@@ -3199,6 +3216,7 @@ c     spline interpolation of rhotor versus rhopol
 
 
       function frhotor_an(rhop)
+      use dierckx, only : splev
       implicit real*8(a-h,o-z)
       parameter(nnr=101,nrest=nnr+4)
       dimension trot(nrest),crot(nrest),rrs(1),ffspl(1)
@@ -4153,6 +4171,7 @@ c
       subroutine equinum_psi(rpsim,zpsim)
       use interp_eqprof, only : rmnm,rmxm,zmnm,zmxm,psiant,psinop,
      . tr,tz,ccspl=>cceq,nsrt,nszt
+      use dierckx, only : fpbisp
 
       implicit real*8 (a-h,o-z)
       parameter(lwrk=8,liwrk=2)
@@ -4263,6 +4282,7 @@ c     here lengths are measured in meters
       end
       
       subroutine sub_derpsi(rpsim,zpsim,nur,nuz,derpsi,cc,lw)
+      use dierckx, only : fpbisp
       use interp_eqprof, only : nr,nz,psia,tr,tz,nsrt,nszt
       
       implicit real*8 (a-h,o-z)
@@ -4288,6 +4308,7 @@ c     here lengths are measured in meters
       end
 
       subroutine equinum_fpol(iderfpol)
+      use dierckx, only : splev,splder
       use interp_eqprof, only : nr,psia,tfp,nsft,cfp,fpolas
       
       implicit real*8 (a-h,o-z)
@@ -4369,7 +4390,7 @@ c
       subroutine psi_raxis(h,r1,r2)
       use interp_eqprof, only : psiant,psinop,zmaxis,nsr=>nsrt,
      . nsz=>nszt,cc=>cceq,tr,tz,nr
-
+      use dierckx, only : profil,sproota
       implicit real*8 (a-h,o-z)
       parameter(mest=4,kspl=3)
       dimension zeroc(mest)
@@ -4421,6 +4442,7 @@ c
       use graydata_anequil, only : dens0,aln1,aln2
       use interp_eqprof, only : psnpp,aad=>aa,bbd=>bb,ccd=>cc,tfn,nsfd,
      . cfn,npp
+      use dierckx, only : splev,splder
 
       implicit real*8 (a-h,o-z)
       dimension xxs(1),ffs(1)
@@ -4478,6 +4500,7 @@ c
       use graydata_anequil, only : te0,dte0,alt1,alt2
       use interp_eqprof, only : psrad,ct,npp
       use utils, only : locate
+      use simplespline, only :spli
 
       implicit real*8 (a-h,o-z)
 c
@@ -4503,6 +4526,7 @@ c
       use graydata_anequil, only : zeffan
       use interp_eqprof, only : psrad,cz,npp
       use utils, only : locate
+      use simplespline, only :spli
 
       implicit real*8 (a-h,o-z)
 c
@@ -5191,6 +5215,7 @@ c
 c
       subroutine valpsispl(rpsi,voli,dervoli,areai,rrii,rbavi,
      .                     bmxi,bmni,fci,intp)
+      use simplespline, only :spli,splid
       use magsurf_data, only : rpstab,cvol,crri,crbav,cbmx,cbmn,
      . carea,cfc,npsi
 
@@ -5222,6 +5247,7 @@ c
 c
       subroutine ratioj(rpsi,ratjai,ratjbi,ratjpli)
       use magsurf_data, only : rpstab,cratja,cratjb,cratjpl,npsi
+      use simplespline, only :spli
 
       implicit real*8 (a-h,o-z)
 
@@ -5909,12 +5935,11 @@ c
 
       subroutine vlambda(alam,psi,fv,dfv)
       use magsurf_data, only : ch,ch01,tjp,tlm,njpt,nlmt,npsi
+      use dierckx, only : fpbisp
 
       implicit real*8 (a-h,o-z)
       parameter (nlam=41)
       parameter(ksp=3,nlest=nlam+ksp+1)
-      
-      external fpbisp
 
       dimension xxs(1),yys(1),ffs(1)
       integer, dimension(:), allocatable :: iwrk

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