scale_current.f90

Go to the documentation of this file.

subroutine scale_current(a, r0, b0, rav, fcirc, kkk)
!----------------------------------------------------------------------
! subroutine to scale the toroidal current profile to match with the
! neoclassical bootstrap current at the edge
! bootstrap current is calculated in subroutine neo using sauter's 
! formula. it is compared with the equilibrium parallel current.
! the input toroidal current is scaled to match the parallel currents.
! the routine can have problems, if the equilibrium parallel current is 
! negative. in that case, it produces flat toroidal current profile in the 
! range where the parallel current is negative. increasing iscalejz can 
! help. 
!
! input parameters (all in mod_dete) : 
!
! iscalejz - the number of scaling loops. 0 = no scaling, normal helena run.
!            default = 0.0
! bsmultip - the amount of the neoclassical bootstrap current 
!            used in the scaling. default = 1.0
! scalemix - how much of the old j_tor is mixed with the new default = 0.0
!----------------------------------------------------------------------

  use itm_constants
  use mod_parameter
  use mod_dat
  use mod_dete
  use mod_gauss
  use mod_mesh
  use mod_nodes
  use mod_profiles
  use mod_helena_io, only : iu6
  use mod_interpolation
  use mod_ne_te_profile
  use helena_spline

  implicit none

  real(r8), intent(in) :: a, r0, b0
  real(r8), dimension(nr), intent(in) :: rav, fcirc
  integer(itm_i4), intent(in) :: kkk

  real(r8), parameter :: zmu0 = 4.e-7_r8 * pi
  real(r8), parameter :: ze = itm_qe
  integer(itm_i4), parameter :: nrbs = 2001

  real(r8), dimension(nrbs) :: psibs
  real(r8), dimension(npts) :: oldzj, oldold, oldoldzj
!    real :: oldoldzjz(npts)
  real(r8), dimension(nr) :: pskn1
  real(r8), dimension(nr) :: p0tmp, ftmp, dptmp, dftmp
  type(spline_coefficients) :: r_spline, c_spline, pp_spline, z_spline, &
   s_spline, ff_spline, q_spline, f_spline
  real(r8), dimension(nr) :: zps, zjj1, zjj2, zjj3, zjj4, zjj5, zjj6, djj5
  real(r8), dimension(nr) :: qq, dq, zjpar
  real(r8), dimension(nr) :: zjbt, zjsig
  real(r8), dimension(3) :: abltg
  real(r8) :: factas
  real(r8) :: radius, cpsurf
  real(r8) :: pscale, rbscale
  real(r8) :: r, s, ws
  real(r8) :: sumq, sumj1, sumj2, sumj3, sumj4, sumj5, sumj6, &
   sumj5r, sumqr, sumb0, sumbop, sumor2
  real(r8) :: x, xr, xs, xrs, xrr, xss
  real(r8) :: y, yr, ys, yrs, yrr, yss
  real(r8) :: ps, psr, pss, psrs, psrr, psss
  real(r8) :: xjac, xjacr, dl, bigr, gradps2, dssdr
  real(r8) :: zjdchi, b02
  real(r8) :: or2av, b02av, gbar
  real(r8) :: ex, ext, denom, dndpsico, dtdpsico, ztmp, defte, tefte
  real(r8) :: conste, consde
  real(r8) :: ss, pwr
  real(r8) :: zni, dni, tte, ti, dtte, dti
  real(r8) :: rr, fcsp, fsp, qsp
  real(r8) :: zznue, zznui, zsigsp, zsigneo, zzjbt, zhjbt
  real(r8) :: zj0, sigz0
!    real :: maxzjpar, maxzbs
  real(r8) :: sumerr, sumerr2
  real(r8) :: dps, fi, pp
  real(r8) :: zj, bsz, sigz, zj2
  real(r8) :: oldzjz, zjzcalc
!  real :: dummyzjz
  real(r8) :: scalejz
  integer(itm_i4) :: iscalest, iscalend
  integer(itm_i4) :: i, j, ngs
  integer(itm_i4) :: typ
  integer(itm_i4) :: n1, n2, n3, n4
  integer(itm_i4) :: ni, nint

!  if (verbosity > 3) write(iu6, *) 'iscalejz = ', iscalejz

  iscalest = scalest * npts
  iscalend = scaleend * npts

  do i = 1, nr - 1
    pskn1(i) = psikn(nr + 1 - i)
  end do 
  pskn1(nr) = 1

  typ = 2
  call allocate_spline_coefficients(r_spline, nr)
  call spline(nr, pskn1, rav, 0._r8, 0._r8, typ, r_spline)
  call allocate_spline_coefficients(c_spline, nr)
  call spline(nr, pskn1, fcirc, 0._r8, 0._r8, typ, c_spline)

  factas = 2._r8
  if (ias == 1) factas = 1._r8

! for r0 = 1.0 and b0 =1.0

  radius = eps
  cpsurf = radius**2 / alfa

  call profiles(p0tmp, ftmp, dptmp, dftmp, a)

  pscale  = eps / alfa**2
  rbscale = 1.0_r8

  do i = 1, nr
    p0tmp(i) = p0tmp(i) * pscale
    dptmp(i) = dptmp(i) * pscale
    ftmp(i) = ftmp(i) * rbscale
    dftmp(i) = dftmp(i) * rbscale
  end do

  r = -1._r8
  do i = 1, nr - 1
    sumq = 0._r8
    sumj1 = 0._r8
    sumj2 = 0._r8
    sumj3 = 0._r8
    sumj4 = 0._r8
    sumj5 = 0._r8
    sumj6 = 0._r8
    sumj5r= 0._r8
    sumqr = 0._r8
    sumb0  = 0._r8
    sumbop = 0._r8
    sumor2 = 0._r8
    do j = 1, np - 1
      n1 = (i - 1) * np + j
      n2 = n1 + 1
      n3 = n2 + np
      n4 = n1 + np
!------------------------------------- 4 point gaussian int. in s -----
      do ngs = 1, 4
        s = xgauss(ngs)
        ws = wgauss(ngs)
        call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
         r, s, x, xr, xs, xrs, xrr, xss)
        call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
         r, s, y, yr, ys, yrs, yrr, yss)
        call interpolation(2, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
         psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps, &
         psr, pss, psrs, psrr, psss)
        xjac = xr * ys - xs * yr
        xjacr = xrr * ys + xr * yrs - xrs * yr - xs * yrr  
        dl = sqrt(xs**2+ys**2)
        bigr = 1._r8 + eps * x
        gradps2 = psr**2 * (xs**2 + ys**2) / xjac**2
        dssdr = abs(psr) / (2._r8 * sqrt(ps))          
!------------------------------------------------ normalisations          
        gradps2 = gradps2 * (cpsurf / radius)**2
        xjac = xjac * radius**2
        bigr = bigr
        dl = radius * dl
        psr = cpsurf * psr
        psrr = cpsurf * psrr
        xjacr = xjacr * radius**2
!---------------------------------------------------------------

        zjdchi = bigr * xjac / abs(psr)

        sumj1 = sumj1 - ws * zjdchi / (gradps2 * bigr**2)
        sumj2 = sumj2 - ws * zjdchi / gradps2
        sumj3 = sumj3 - ws * zjdchi * bigr**2 / gradps2
        sumj4 = sumj4 - ws * zjdchi / bigr**2
        sumj5 = sumj5 - ws * zjdchi
        sumj6 = sumj6 - ws * zjdchi * gradps2 / bigr**2 
        sumq  = sumq - ws * xjac / (bigr * abs(psr))

        sumqr = sumqr + psrr * xjac / ((psr**2) * bigr) * ws
        sumqr = sumqr - xjacr / (bigr * psr) * ws
        sumqr = sumqr + xjac * eps * xr / ((bigr**2) * psr) * ws

        sumj5r = sumj5r + xjac  * bigr * psrr / (psr**2) * ws
        sumj5r = sumj5r - xjacr * bigr / psr * ws
        sumj5r = sumj5r - xjac * eps *xr / psr * ws

!----------------------------------------- d_nc from Hegna's

        b02 = (ftmp(nr - i + 1) / bigr)**2 + gradps2 / bigr**2

        sumb0  = sumb0 - ws * zjdchi * b02
        sumbop = sumbop - ws * zjdchi * b02 / gradps2
        sumor2 = sumor2 - ws * zjdchi / bigr**2

      end do
    end do
    ni = nr - i + 1
    zps(ni) = sqrt(ps)
    zjj1(ni) = factas * sumj1 / (2._r8 * pi)
    zjj2(ni) = factas * sumj2 / (2._r8 * pi)
    zjj3(ni) = factas * sumj3 / (2._r8 * pi)
    zjj4(ni) = factas * sumj4 / (2._r8 * pi)
    zjj5(ni) = factas * sumj5 / (2._r8 * pi)
    zjj6(ni) = factas * sumj6 / (2._r8 * pi)

    djj5(ni) = factas * sumj5r / dssdr / (2._r8 * pi)

    qq(ni) = factas * ftmp(ni) * sumq / (2._r8*pi)
    dq(ni) = dftmp(ni) * sumq + ftmp(ni) * sumqr / dssdr
    dq(ni) = dq(ni) * factas / (2._r8 * pi)

    or2av = factas * sumor2 / (2._r8 * pi) / zjj5(ni)
    b02av = factas * sumb0 / (2._r8 * pi) / zjj5(ni)
    gbar = factas * sumbop / (2._r8 * pi) 

    zjpar(ni) = (-dptmp(ni) * ftmp(ni) - dftmp(ni) * b02av) &
     / ( 2._r8 * cpsurf * sqrt(ps))
    zjpar(ni) = zjpar(ni) / (zmu0 * r0 / b0**2)
!    if (verbosity > 3) write(iu6, *) ni, qq(ni)
  end do
  qq(nr) = 2._r8 * qq(nr - 1) - qq(nr - 2)
  qq(1) = 2._r8 * qq(2) - qq(3)

  radius = eps * r0
  cpsurf = radius**2 * b0 / alfa

!  if (verbosity > 3) write(iu6, *) r0, b0, alfa, cpsurf

  call profiles(p0tmp, ftmp, dptmp, dftmp, a)

  pscale  = b0**2 * eps / (zmu0 * alfa**2)
  rbscale = b0 * r0

  do i = 1, nr
    p0tmp(i) = p0tmp(i) * pscale
!TODO: This was disabled in revision 161. scale_current has to be re-written
!    if (p%type == 15) then 
!      p0tmp(i) = p0tmp(i) + p%coeff(2) * 1000._r8
!    end if
    dptmp(i) = dptmp(i) * pscale
    ftmp(i) = ftmp(i) * rbscale
    dftmp(i) = dftmp(i) * rbscale
  end do

  call allocate_spline_coefficients(pp_spline, nr)
  call spline(nr, pskn1, dptmp, 0._r8, 0._r8, typ, pp_spline)

  if (te%shape > 0) then
    call ne_te_profile(ne, ne%psi_0, defte, dndpsico, consde)
    call ne_te_profile(te, te%psi_0, tefte, dtdpsico, conste)
  end if

  call allocate_spline_coefficients(z_spline, nr)
  call allocate_spline_coefficients(s_spline, nr)
  call allocate_spline_coefficients(ff_spline, nr)
  call allocate_spline_coefficients(q_spline, nr)
!------------------------------------ to be elaborated
  if (te%shape == 0) then
    do i = 2, nr - 1
      ps = psi(4 * (nr - i) * np + 1)
      ss = sqrt(ps)

      pwr = 1._r8 / (etaei + 1._r8)

      zni = zn0 * (p0tmp(i) / p0tmp(1))**pwr
      dni = zn0 * pwr * abs((p0tmp(i) / p0tmp(1)))**(pwr - 1._r8) &
       * dptmp(i) / p0tmp(1) / (2._r8 * ss)
      tte = rpe * p0tmp(i) / (ze * zni * 1e19_r8)
      ti = (1._r8 - rpe) * p0tmp(i) / (ze * zni * 1e19_r8)

      dtte = rpe / (ze * 1e19_r8) * (dptmp(i) / (2._r8 * ss * zni) &
       - p0tmp(i) * dni / zni**2)
      dti = (1._r8 - rpe) / (ze * 1e19_r8) * (dptmp(i) / (2._r8 &
       * ss * zni) - p0tmp(i) * dni / zni**2)
      dtte = dtte / cpsurf
      dti = dti / cpsurf
      dni = dni / cpsurf
      rr = r0 * rav(i)

      call neo_helena(fcirc(i), tte, ti, zni, zeff, ftmp(i), qq(i), rr, &
       eps * ss, dni, dtte, dti, zznue, zznui, zsigsp, zsigneo, &
       zzjbt, zhjbt)
      zjbt(i) = zzjbt
      if (neosig) then
        zjsig(i) = zsigneo
      else
        zjsig(i) = zsigsp
      end if
      if (i > 2) zjsig(i) = zjsig(i) / zjsig(2)
    end do
    stop
    zjsig(2) = 1.0_r8
    zjsig(1) = 2._r8 * zjsig(2) - zjsig(3)
    zjsig(nr) = 2._r8 * zjsig(nr - 1) - zjsig(nr - 2)
    zjbt(1) = zjbt(2)
    zjbt(nr) = zjbt(nr - 1)
    call spline(nr, pskn1, zjbt, 0._r8, 0._r8, typ, z_spline)
    call spline(nr, pskn1, zjsig, 0._r8, 0._r8, typ, s_spline)
  else
    call spline(nr, pskn1, ftmp, 0._r8, 0._r8, typ, ff_spline)
    call spline(nr, pskn1, qq, 0._r8, 0._r8, typ, q_spline)
    do i = 2, nrbs - 1
      psibs(i) = dble(i - 1) / dble(nrbs)
      ps = psibs(i)
      ss = sqrt(ps)
      if (ps >= ne%psi_0) then
        call ne_te_profile(ne, ps, zni, dni)
      else
        call ne_te_profile(ne, ps, zni, dni, consde, dndpsico)
!TODO: factor fde missing in zni (same in si_output) !!!
      end if
      dni = -dni
!     te given in ev (profiles defined in kev)
      if (ps >= te%psi_0) then
        call ne_te_profile(te, ps, tte, dtte)
      else
        call ne_te_profile(te, ps, tte, dtte, conste, dtdpsico)
      end if              
      dtte = -dtte * 1000.0_r8
      ti = tte
      dti = dtte
      dtte = dtte / cpsurf
      dti = dti / cpsurf
      dni = dni / cpsurf

!      if (verbosity > 3) write(iu6, *) ps, te, zni, dte, dni, conste, consde
      rr = r0 * spwert(nr, ps, r_spline, pskn1, abltg, 0)
      fcsp = spwert(nr, ps, c_spline, pskn1, abltg, 0)
      fsp = spwert(nr, ps, ff_spline, pskn1, abltg, 0)
      qsp = spwert(nr, ps, q_spline, pskn1, abltg, 0)
!      if (i < 200) write(iu6, *) ps, fcsp, te, ti, zni, qsp, dni, dte, dti

      call neo_helena(fcsp, tte, ti, zni, zeff, fsp, qsp, rr, eps * ss, dni, &
       dtte, dti, zznue, zznui, zsigsp, zsigneo, zzjbt, zhjbt)
!      if (i == nrbs - 3) then
!       if (verbosity > 3) &
!        write(iu6, *) i, ps, fcirc(i), te, ti, zni, qsp, dni, dte, dti, &
!        zzjbt, zznue, conste, dtdpsico               
!      end if
      zjbt(i) = zzjbt
      if (neosig) then
        zjsig(i) = zsigneo
      else
        zjsig(i) = zsigsp
      end if
      if (i > 2) zjsig(i) = zjsig(i) / zjsig(2)
!      if (i < 200) write(iu6, *) ps, te, ti, zni, dni, dte, dti, &
!       zsigsp, zzjbt
    end do  
    zjsig(2) = 1.0_r8
    zjsig(1) = 2._r8 * zjsig(2) - zjsig(3)
    zjsig(nrbs) = 2._r8 * zjsig(nrbs - 1) - zjsig(nrbs - 2)
!    zjbt(1) = zjbt(2)
    zjbt(nrbs) = 2._r8 * zjbt(nrbs - 1) - zjbt(nrbs - 2)
    zjbt(1) = 0.0_r8
    psibs(1) = 0.0_r8
    psibs(nrbs) = 1.0_r8
    call spline(nrbs, psibs, zjbt, 0._r8, 0._r8, typ, z_spline)
    call spline(nrbs, psibs, zjsig, 0._r8, 0._r8, typ, s_spline)
!    do i = 1, nrbs
!      if (verbosity > 3) write(iu6, *) psibs(i), zjbt(i), zjsig(i)
!    end do
!    stop
  end if
  zjpar(1) = 2._r8 * zjpar(2) - zjpar(3)
!  zjpar(nr) = 2 * zjpar(nr - 1) - zjpar(nr - 2)
  zj0 = zjpar(1)
  sigz0 = zjsig(1)
  call allocate_spline_coefficients(f_spline, nr)
  call spline(nr, pskn1, zjpar, 0._r8, 0._r8, typ, f_spline)
!  if (verbosity > 3) write(iu6, *) 'iscalejz = ', iscalejz
  if (iscalejz > 0) then
    if (verbosity > 3) write(iu6, *) 'scaling current'
    if (verbosity > 3) write(iu6, *) 'scalemix = ', scalemix
!    maxzjpar = 0._r8
!    maxzbs = 0._r8
!    do i = scaleend * nr, nr - 1
!      if (zjpar(i) > maxzjpar) then
!        maxzjpar = zjpar(i)
!      end if
!      if (zjbt(i) > maxzbs) then
!        maxzbs = zjbt(i)
!      end if
!      if (verbosity > 3) write(iu6, *) i, zjpar(i)
!    end do
!    if (verbosity > 3) write(iu6, *) 'scalejz = ', scalejz, maxzbs, maxzjpar
!    if (verbosity > 3) write(iu6, *) iscalest, iscalend
    sumerr = 0.0_r8
    sumerr2 = 0.0_r8
    if (verbosity > 3) &
     write(iu6, *) "flux   j_tor(i)_cal      j_tor(i)      oldzjz      ", &
     "    scalejz          bsz         sigz         zj          pp"
    dps = 1._r8 / dble(npts - 1)
    do i = 1, npts
      fi = dble(i - 1) / dble(npts - 1)
      if (i == 1) then
        zj = zjpar(1)
        bsz = zjbt(1)
        sigz = zjsig(1)
      else
        zj = spwert(nr, fi, f_spline, pskn1, abltg, 0)
        nint = int((npts - 1) * fi) + 1
        zj2 = zjpar(nint) + (fi - dps * (nint - 1)) / dps &
         * (zjpar(nint + 1) - zjpar(nint))
!        if (verbosity > 3) write(iu6, *) nint, fi, zj, zj2
        if (te%shape == 0) then
          bsz = spwert(nr, fi, z_spline, pskn1, abltg, 0)
          sigz = spwert(nr, fi, s_spline, pskn1, abltg, 0)
        else
          bsz = spwert(nrbs, fi, z_spline, psibs, abltg, 0)
          sigz = spwert(nrbs, fi, s_spline, psibs, abltg, 0)
        end if             
      end if

      pp = -spwert(nr, fi, pp_spline, pskn1, abltg, 0)
!      if (fi > 0.97_r8) then
       bsz = bsz * bsmultip
!      end if
!      end do
!      do i = iscalend, npts - 1
!        fi = dble(i - 1) / dble(npts - 1)
!        zj = spwert(nr, fi, f_spline, pskn1, abltg, 0)
!        if (te%type == 0) then
!          bsz = spwert(nr, fi, z_spline, pskn1, abltg, 0)
!          sigz = spwert(nr, fi, s_spline, pskn1, abltg, 0)
!        else
!          bsz = spwert(nrbs, fi, z_spline, psibs, abltg, 0)
!          sigz = spwert(nrbs, fi, s_spline, psibs, abltg, 0)
!        end if 
!        bsz = bsz * bsmultip
      oldzjz = j_tor(i)
!        if (i >= iscalend) then
      scalejz = bsz / zj
!        sumerr = sumerr + abs(scalejz * (1._r8 - zj / zj0 * sigz &
!         / sigz0) - 1._r8)
      j_tor(i) = sigz + j_tor(i) * scalejz ! *(1-zj/zj0*sigz/sigz0)
      zjzcalc = j_tor(i)
!        else
!          j_tor(i) = sigz
!        end if
!        if (j_tor(i) < 0._r8) j_tor(i) = j_tor(i - 1)
!        if (kkk > 2) then 
!          dummyzjz = ((oldzjz - oldold(i) * (sigz * zj + bsz &
!           - oldzj(i))) / (oldzj(i) - oldoldzj(i))) + oldzjz
!        end if
!        if (j_tor(i) < 0._r8) j_tor(i) = oldzjz * 2._r8
!        if ((zj - bsz) / zj0 < sigz .and. zj > bsz) then
!          j_tor(i) = oldzjz * sigz / (zj - bsz) * zj0
!        end if
      if (zj < 0._r8) then
        if (j_tor(i - 1) > j_tor(i) .and. j_tor(i - 1) > oldzjz) then
          j_tor(i) = j_tor(i - 1)
!          if (verbosity > 3) &
!           write(iu6, *) 'j_tor(',i,') < 0 !!', j_tor(i), j_tor(i-1), oldzjz
        else
          j_tor(i) = oldzjz * 2.0_r8
        end if
      end if
      zjzcalc = j_tor(i)
!      if ((zj - bsz) / zj0 > sigz) then
!        if (verbosity > 3) write(iu6, *) 'j_tor larger than sigz'
!        j_tor(i) = oldzjz * sigz / (zj - bsz) * zj0
!      endif
      if (i > 3 * npts / 4) then
        sumerr = sumerr + abs(oldzjz / j_tor(i) - 1.0_r8)
        sumerr2 = sumerr2 + abs(bsz / (zj - sigz / sigz0 * zj0) &
         - 1.0_r8)
      end if
      if (i > 1) then
!        if (((oldzjz .lt. oldold(i - 1) .and. oldzjz .lt. oldold(i &
!         + 1) .and. j_tor(i) .gt. oldzjz) .or. (oldzjz .gt. oldold(i &
!         - 1). and. oldzjz .gt. oldold(i + 1)) .and. j_tor(i) &
!         .lt. oldzjz) .and. i .lt. npts) then
!          j_tor(i) = scalemix / 2 * oldzjz + (1.0 - scalemix / 2) &
!           * j_tor(i) / j_tor(1)
!        else
        j_tor(i) = scalemix * oldzjz + (1.0_r8 - scalemix) * j_tor(i) &
         / j_tor(1)
!       endif
      end if
      if (i > iscalend .or. i == 1) then            
        if (verbosity > 3) &
         write(iu6, 10) fi, zjzcalc, j_tor(i), oldzjz, scalejz, bsz, &
         sigz, zj, pp
      end if
      oldold(i) = oldzjz
      oldoldzj(i) = oldzj(i)
      oldzj(i) = zj + sigz
    end do
    j_tor(1) = 1._r8
    if (verbosity > 3) &
     write(iu6, *) 'average error in current profile vs. ', &
     'bootstrap current: ', sumerr2 / (npts / 4)
  
    if (verbosity > 3) write(iu6, *) 'average error in current iteration: ', &
     sumerr / (npts / 4)
!    j_tor(npts) = 2 * j_tor(npts - 1) - j_tor(npts - 2)
!    do i = iscalest, iscalend - 1
!      oldzjz = j_tor(i)
!      j_tor(i) = j_tor(iscalest) + (j_tor(iscalend) - j_tor(iscalest)) &
!       * (i - iscalest) / (scaleend - scalest) / dble(npts)
!      if (verbosity > 3) write(iu6, *) i, j_tor(i), oldzjz            
!    end do
  end if
!  if ((kkk > 0.5_r8 * iscalejz) .and. (iscalejz >= 10) &
!   .and. (scalemix > 0.5_r8)) then
!    scalemix = 0.5_r8
!  end if
!  if ((kkk > 0.8_r8 * iscalejz) .and. (iscalejz >= 10) &
!   .and. (scalemix > 0.2_r8)) then
!    scalemix = 0.2_r8
!  end if

   call deallocate_spline_coefficients(r_spline)
   call deallocate_spline_coefficients(c_spline)
   call deallocate_spline_coefficients(ff_spline)
   call deallocate_spline_coefficients(f_spline)
   call deallocate_spline_coefficients(z_spline)
   call deallocate_spline_coefficients(q_spline)
   call deallocate_spline_coefficients(s_spline)
   call deallocate_spline_coefficients(pp_spline)

 10   format(f6.4, 9es14.6)
  
  return
end subroutine scale_current