!-----------------------------------------------------------------------------
subroutine kdsaw_pre_init_r4(zfp_r4,imix_opt)
  !  mixing model select (imix_opt)
  !    1 -- standard Kadomtsev
  !    2 -- Kadomtsev with guaranteed flat n & T profiles post-sawtooth
  !    3 -- standard Porcelli
  !    4 -- Porcelli with guaranteed flat n & T profiles post-sawtooth
  !
  !  R4 porcelli mixing fraction preset (zfp_r4)
  !    1.0 --> full Kadomtsev-sized q=1 mixing region (all the way to x=0).
  !    0.9 --> q=1 island is 90% of Kadomtsev island(in x flux coordinate).
  !    0.0 --> zero width q=1 island.
  !  Even with the Porcelli model, the axial region inside the q=1 island
  !  gets  mixed.

  implicit NONE
  real, intent(in) :: zfp_r4
  integer, intent(in) :: imix_opt

  real*8 :: zfp

  zfp = zfp_r4
  call kdsaw_pre_init(zfp,imix_opt)

end subroutine kdsaw_pre_init_r4

!-----------------------------------------------------------------------------
subroutine kdsaw_pre_init(zfp,imix_opt)
  !  mixing model select (imix_opt)
  !    1 -- standard Kadomtsev
  !    2 -- Kadomtsev with guaranteed flat n & T profiles post-sawtooth
  !    3 -- standard Porcelli
  !    4 -- Porcelli with guaranteed flat n & T profiles post-sawtooth
  !
  !  porcelli mixing fraction preset (zfp)
  !    1.0 --> full Kadomtsev-sized q=1 mixing region (all the way to x=0).
  !    0.9 --> q=1 island is 90% of Kadomtsev island(in x flux coordinate).
  !    0.0 --> zero width q=1 island.
  !  Even with the Porcelli model, the axial region inside the q=1 island
  !  gets  mixed.

  use kdsaw_module

  implicit NONE
  real(fp), intent(in) :: zfp
  integer, intent(in) :: imix_opt

  f_porcelli = max(ZERO,min(ONE,zfp))
  mix_method = max(1,min(4,imix_opt))

end subroutine kdsaw_pre_init

!-----------------------------------------------------------------------------
subroutine kdsaw_dxi_qjump_r4(zdxi)

  !  set finite thickness of current sheets (this will have an upper bound of
  !  1/20 the sawtooth mixing region width).

  real, intent(in) :: zdxi
  real*8 :: zdxi_r8

  zdxi_r8 = zdxi
  call kdsaw_dxi_qjump(zdxi_r8)

end subroutine kdsaw_dxi_qjump_r4

!-----------------------------------------------------------------------------
subroutine kdsaw_dxi_qjump(zdxi)

  !  set finite thickness of current sheets (this will have an upper bound of
  !  1/20 the sawtooth mixing region width).

  use kdsaw_module
  implicit NONE

  real(fp), intent(in) :: zdxi

  real(fp), parameter :: zmax_dxi = 0.1_fp

  dxi_qjump = min(zmax_dxi,zdxi)
  if(dxi_qjump.lt.zero) dxi_qjump=zero

end subroutine kdsaw_dxi_qjump

!-----------------------------------------------------------------------------
subroutine kdsaw_dxi_qjump_get_r4(zdxi)

  !  fetch finite thickness of current sheets

  real, intent(out) :: zdxi
  real*8 :: zdxi_r8

  call kdsaw_dxi_qjump_get(zdxi_r8)
  zdxi = zdxi_r8

end subroutine kdsaw_dxi_qjump_get_r4

!-----------------------------------------------------------------------------
subroutine kdsaw_dxi_qjump_get(zdxi)

  !  fetch finite thickness of current sheets

  use kdsaw_module
  implicit NONE

  real(fp), intent(out) :: zdxi

  zdxi = dxi_qjump

end subroutine kdsaw_dxi_qjump_get

!-----------------------------------------------------------------------------
subroutine kdsaw_porcelli_get_r4(zfp_r4,imix_opt)

  !  get porcelli mixing fraction

  implicit NONE
  real, intent(out) :: zfp_r4
  integer, intent(out) :: imix_opt

  real*8 :: zfp

  call kdsaw_porcelli_get(zfp,imix_opt)
  zfp_r4 = zfp 

end subroutine kdsaw_porcelli_get_r4

!-----------------------------------------------------------------------------
subroutine kdsaw_porcelli_get(zfp,imix_opt)

  !  get porcelli mixing fraction 

  use kdsaw_module

  implicit NONE
  real(fp), intent(out) :: zfp
  integer, intent(out) :: imix_opt

  zfp = f_porcelli
  imix_opt = mix_method

end subroutine kdsaw_porcelli_get

!-----------------------------------------------------------------------------
subroutine kdsaw_init_r4(ilun_msg,izones,zxi_r4,ziotb_r4,zdvol_r4,iimix,ierr)
  !  real interface... see kdsaw_init (below)
  !  dmc 24 Jan 2002

  implicit NONE

  integer, intent(in) :: ilun_msg,izones
  real, intent(in) :: zxi_r4(izones+1),ziotb_r4(izones+1),zdvol_r4(izones)
  integer, intent(out) :: iimix, ierr

  !-------------------

  integer :: ilun
  real*8 :: zxi(izones+1),ziotb(izones+1),zdvol(izones)

  !-------------------

  if(izones.le.5) then
     ilun=abs(ilun_msg)
     write(ilun,1001) izones
 1001 format(' ?kdsaw_init_r4: izones=',i10,' invalid, require izones.ge.5')
     ierr=1
     return
  endif

  zxi = zxi_r4
  ziotb = ziotb_r4
  zdvol = zdvol_r4

  call kdsaw_init(ilun_msg,izones,zxi,ziotb,zdvol,iimix,ierr)

  return
end subroutine kdsaw_init_r4

!-----------------------------------------------------------------------------
subroutine kdsaw_init(ilun_msg,izones,zxi,ziotb,zdvol,iimix,ierr)

  !  initialize Kadomtsev sawtooth mixing model (real*8 interface)
  !  dmc 24 Jan 2002

  use kdsaw_module
#ifdef __DEBUG
  use ezcdf
#endif

  implicit NONE

  !  input data:

  integer, intent(in) :: ilun_msg
                         ! fortran i/o channel for error messages

  integer, intent(in) :: izones
                         ! number of plasma zones
  real(fp), dimension(izones+1) :: zxi
                         ! zone *bdy* locations in normalized sqrt(tor.flux)
                         ! in strict ascending order with...
                         ! zxi(1)=0, zxi(izones+1)=1 expected
                         ! zxi(j) is inner bdy, zxi(j+1) outer bdy of j'th zone

  real(fp), dimension(izones+1) :: ziotb
                         ! iota(bar) (1/q) profile at zone bdys

  real(fp), dimension(izones) :: zdvol
                         ! "volume" of the zones 
                         ! in arbitrary units, only the relative values
                         ! amongst zones are important

  ! output data:

  integer, intent(out) :: iimix
                         ! iimix gives index to last zone involved in
                         ! Kadomtsev mixing, or 0 if none (not an error).

  integer, intent(out) :: ierr
                         ! error flag-- error in input data,
                         !   such as, izones.lt.5, or, zxi out of order.

  !---------------------------------------------------------------

  integer i,icdf

  !---------------------------------------------------------------

#ifdef __DEBUG
  if(ilun_msg.ge.0) then
     call cdf_open(icdf,'kdsaw_debug.cdf','w')
     call cdf_define(icdf,'nzones',izones)
     call cdf_define(icdf,'xi',zxi)
     call cdf_define(icdf,'xiotb',ziotb)
     call cdf_define(icdf,'dvol',zdvol)
     call cdf_write(icdf,'nzones',izones)
     call cdf_write(icdf,'xi',zxi)
     call cdf_write(icdf,'xiotb',ziotb)
     call cdf_write(icdf,'dvol',zdvol)
     call cdf_close(icdf)
  endif
#endif
  if(kdsaw_flag) then
     deallocate(xi,xiotbp,xiotb,helflx,dvol)
     deallocate(iord,izonec,zxipre,inzmix,izmix,zzmix)
     deallocate(inzmxp,izmixp,zzmixp)
     kdsaw_flag=.FALSE.
  endif

  iimix=0
  ierr=0

  !  error checks 

  if(izones.le.5) then
     write(abs(ilun_msg),1001) izones
     ierr=1
     return
  endif
1001 format(' ?kdsaw_init: izones=',i10,' invalid, require izones.ge.5')

  do i=1,izones
     if(zxi(i).ge.zxi(i+1)) then
        write(abs(ilun_msg),1002) i,zxi(i),zxi(i+1)
        ierr=2
        return
     endif
  enddo
1002 format(' ?kdsaw_init: in zone i=',i5,' zxi outer bdy < inner bdy:'/ &
          ' zxi(i)=',1pd11.4,' zxi(i+1)=',1pd11.4)

  if(abs(zxi(1)).gt.1.0e-5_fp) then
     write(abs(ilun_msg),1003) zxi(1)
     ierr=3
     return
  endif
  if(abs(zxi(izones+1)-1).gt.1.0d-5) then
     write(abs(ilun_msg),1004) zxi(izones+1)
     ierr=4
     return
  endif
1003 format(' ?kdsaw_init: inner bdy must be zero, but zxi(1)=',1pd11.4)
1004 format(' ?kdsaw_init: outer bdy must be one, but zxi(izones+1)=',1pd11.4)

  do i=1,izones
     if(ziotb(i)*ziotb(izones+1).le.zero) then
        write(abs(ilun_msg),1005) i,ziotb(i),ziotb(izones+1)
        ierr=5
        return
     endif
     if(zdvol(i).le.0) then
        write(abs(ilun_msg),1006) i,zdvol(i)
        ierr=6
        return
     endif
  enddo
1005 format(' ?kdsaw_init: iota(bar) profile sign changes:'/ &
          '  ziotb(',i5,')=',1pd11.4,'; ziotb(izones+1)=',1pd11.4)
1006 format(' ?kdsaw_init: zero or negative volume element:'/ &
          '  zdvol(',i5,')=',1pd11.4)

!-------------------------------------------------------
!  OK to set up module...
!
  kdsaw_flag = .TRUE.

  nonlin=ilun_msg

  nzones=izones
  mj=nzones+3
  lcentr=2
  lcp1=lcentr+1
  ledge=lcentr+nzones-1
  lep1=ledge+1

  !  old TRANSP-style array allocation & setup

  allocate(xi(mj,2)); xi=zero

  allocate(xiotb(mj,2),xiotbp(mj,2),helflx(mj))
  xiotb=zero; xiotbp=zero; helflx=zero

  allocate(dvol(mj,2)); dvol=zero
  
  !  grid

  xi(lcentr:lep1,1)=zxi
  xi(lcentr,1)=zero
  xi(lep1,1)=one
  xi(lcentr:ledge,2)=(xi(lcentr:ledge,1)+xi(lcp1:lep1,1))/two

  !  pre-sawtooth iota(bar)

  xiotbp(lcentr:lep1,1)=ziotb
  xiotbp(lcentr:ledge,2)=(xiotbp(lcentr:ledge,1)+xiotbp(lcp1:lep1,1))/two

  xiotb=xiotbp          ! post-sawtooth iota(bar) = pre-sawtooth iota(bar)
                        ! initially...

  tflux=one             ! toroidal flux is "normalized".

  !  zone volumes (can be any units)

  dvol(lcentr:ledge,2)=zdvol
  dvol(lep1,2)=two*dvol(ledge,2)-dvol(ledge-1,2)
  dvol(lcp1:lep1,1)=(dvol(lcentr:ledge,2)+dvol(lcp1:lep1,2))/two

  !  compute helical flux = poloidal flux - [polodal flux if q=1 everywhere]

  helflx(lcentr)=zero
  do i=lcentr,ledge
     helflx(i+1) = helflx(i) + &
          tflux*(xiotbp(i,2)-one)*xi(i,2)*(xi(i+1,1)-xi(i,1))/xpi
  enddo

  !-----------------------------------------
  !  allocate "SCRATCH" sawtooth model space

  imix=0
  imax=0

  allocate(iord(mj)); iord=0
  allocate(izonec(mj)); izonec=0

  allocate(zxipre(2,mj)); zxipre=zero

  allocate(izmix(mj,mj)); izmix=0
  allocate(inzmix(mj)); inzmix=0

  allocate(zzmix(mj,mj)); zzmix=zero

  !-----------------------------------------
  !  allocate "SCRATCH" sawtooth model space

  allocate(inzmxp(mj)); inzmxp=0
  allocate(izmixp(mj,mj)); izmixp=0
  allocate(zzmixp(mj,mj)); zzmixp=zero

  !-----------------------------------------
  !  OK... is there a sawtooth?

  call kdsaw_hck(ierr)

  iimix=max(0,imix-lcentr+1)
  if(iimix.eq.0) then
     return      ! no q=1 surface; no sawtooth (or ierr was set in kdsaw_hck).
  endif

  !  There is a sawtooth... derive mixing information from
  !  Helical Flux profile (mix_method.eq.1 or 2) only...

  if(mix_method.le.2) then
     call kdsaw_hdr(ierr)
     if(ierr.ne.0) then
        write(abs(nonlin),*) ' ?kdsaw_init: algorithm error in kdsaw_hdr.for!'
        ierr=99
        return
     endif
  endif

  !  Compute mixing factors
  !  i.e. which presawtooth zones mix in to form each post-sawtooth zone

  call kdsaw_mix(ierr)
  if(ierr.ne.0) then
     write(abs(nonlin),*) ' ?kdsaw_init: algorithm error in kdsaw_mix.for!'
     ierr=99
     return
  endif

  !  Compute distribution probability factors
  !  i.e. probability distribution of particles from each pre-sawtooth zone
  !  being mapped to the various post-sawtooth zones

  call kdsaw_inv(ierr)
  if(ierr.ne.0) then
     write(abs(nonlin),*) ' ?kdsaw_init: algorithm error in kdsaw_inv.for!'
     ierr=99
     return
  endif

  !  All done...

  return
end subroutine kdsaw_init
