real(fp)function zf1d(zzh,zf0,zf1,zf2,zf3)
  use iso_c_binding, only: fp => c_double
  real(fp) :: zzh,zf0,zf1,zf2,zf3
  zf1d = ((zf3*zzh + zf2)*zzh + zf1)*zzh + zf0
end function zf1d

real(fp) function zdf1d(zzh,zf1,zf2,zf3)
  use iso_c_binding, only: fp => c_double
  real(fp) :: zzh,zf1,zf2,zf3
  zdf1d = (3.0_fp*zf3*zzh + 2.0_fp*zf2)*zzh + zf1      ! df/xi
end function zdf1d

real(fp) function zd2f1d(zzh,zf2,zf3)
  use iso_c_binding, only: fp => c_double
  real(fp) :: zzh,zf2,zf3
  zd2f1d = 6.0_fp*zf3*zzh + 2.0_fp*zf2                  ! d2f/dxi2
end function zd2f1d

real(fp) function zpsint(zzprev,zzkonst,zzx0,zzh,zf0,zf1,zf2,zf3) 
  use iso_c_binding, only: fp => c_double
  real(fp) :: zzh,zf0,zf1,zf2,zf3,zzprev,zzkonst,zzx0
  !
  !  poloidal flux integral formula, from iota(bar) spline.
  !  based on d(psi)/d(xi) = iota(bar)*xi*(tor.flux.at.bdy)/xpi
  !
  zpsint = zzprev   + &
       zzkonst * zzh * (zzx0*zf0 + zzh*((zf0 + zzx0*zf1)/2.0_fp+ &
       zzh*((zf1+zzx0*zf2)/3.0_fp+ zzh*((zf2 + zzx0*zf3)/4.0_fp+ &
       zzh*zf3/5.0_fp))))
 
end function zpsint

real(fp) function  zlope0(delta,f0,f0pj,f0p2j,f0p3j)
  !    quadratic extrapolation of slope at xi=0.0 (cubic derivative assumed
  !    constant near edges---second derivative is continuous as in splines)
  !    extrapolated from first interior point (f0pj)
  !    first and second derivatives evaluated at first interior point f0pj
  !    cubic derivative evaluated halfway between f0pj and f0p2j and
  !    assumed to be the same at f0pj
  !
  use iso_c_binding, only: fp => c_double
  real(fp) :: delta,f0,f0pj,f0p2j,f0p3j
  zlope0 = &
       (f0p2j-f0)/(2.0_fp*delta) -&	!constant term &
       (f0+f0p2j-2.0_fp*f0pj)/delta +&	!linear term &
       0.5_fp*(3.0_fp*f0pj+f0p3j-f0-3.0_fp*f0p2j)/delta !quadratic term
end function zlope0

real(fp) function  zlope1(delta,f1,f1mj,f1m2j,f1m3j)
  !    quadratic extrapolation of slope at xi=zdx*(ixi-1)
  use iso_c_binding, only: fp => c_double
  real(fp) :: delta,f1,f1mj,f1m2j,f1m3j
  zlope1 = &
       (f1-f1m2j)/(2.0_fp*delta) + &	!constant term &
       (f1+f1m2j-2.0_fp*f1mj)/delta + &	!linear term &
       0.5_fp*(3.0_fp*f1m2j+f1-f1m3j-3.0_fp*f1mj)/delta	!quadtaric term
end function zlope1
  !
