module kdsaw_module
  
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: zero, half, one, two
  use physconst_mod, only : xpi=>pi 
  IMPLICIT NONE
  SAVE

  ! Kadomtsev sawtooth module -- internal data
  !
  ! dmc 23 Jan 2002 -- this is the Kadomtsev sawtooth code from TRANSP
  ! the module contains local copies of the portions of TRCOM which
  ! it previously used, and, its own internal COMMON (formerly called
  ! "SCRATCH").
  !

  logical :: kdsaw_flag = .FALSE.
                           ! flag that module has been initialized

  !---------------------------------------
  ! Porcelli mixing fraction

  real(fp) :: f_porcelli = one  ! 1 means full Kadomtsev, <1 means Porcelli
  integer :: mix_method = 1     ! 1 means standard, 2 means ergodic everywhere
  real(fp) :: dxi_qjump = one   ! specify finite (in xi coordinate) current
                                ! sheet thickness.  "qjump" refers to a
                                ! discontinuity in q or 1/iota(bar) that
                                ! is equivalent to a current sheet.

  ! standard mixing involves zone matching by Helical Flux (Kadomtsev) or
  ! two regions of ergodic mixing (Porcelli), the two regions being the
  ! island around the q=1 surface and the axial region enclosed by the
  ! island; mix_method=2 means ergodic mixing over the entire mixing region.
  ! Ergodic mixing redistributes particles randomly with probability 
  ! proportional to the toroidal flux subtended by each zone in the mixing 
  ! region.
  !---------------------------------------

  ! originally from TRANSP's TRCOM...

  integer :: nzones = 0    ! number of active zones
  integer :: mj            ! dimension of radial arrays
  integer :: lcentr        ! index to innermost zone
  integer :: lcp1          ! index to outer bdy of innermost zone, lcentr+1
  integer :: ledge         ! index to outermost zone
  integer :: lep1          ! index to outer bdy of outermost zone, ledge+1

  integer :: nonlin = 6    ! i/o channel for messages

  real(fp), dimension(:,:), allocatable :: xi
                           ! grid (normalized sqrt toroidal flux)

  real(fp) :: tflux          ! enclosed toroidal flux, Wb

  real(fp), dimension(:,:), allocatable :: xiotbp
                           ! pre-sawtooth iota(bar) profile

  real(fp), dimension(:), allocatable :: helflx
                           ! pre-sawtooth helical flux profile

  real(fp), dimension(:,:), allocatable :: xiotb
                           ! Kadomtsev idel post-sawtooth iota(bar) profile

  real(fp), dimension(:,:), allocatable :: dvol
                           ! zone volumes (to convert densities)

  real(fp) :: xswmix         ! mixing xi-radius
  real(fp) :: xswinv         ! q=1 (inversion) xi-radius
  real(fp) :: xswmin         ! inside extent of q=1 island (0 for Kadomtsev)

  real(fp) :: zvolq1         ! volume of q=1 island
  real(fp) :: zvolax         ! volume of axial region inside q=1 island
  real(fp) :: zvolmix        ! total volume of mixing region, zvolq1+zvolax.

  real(fp) :: zfrmix         ! fraction of volume of zone (imix) involved
  real(fp) :: zfrmin         ! fraction of volume of zone (imin) in axial region

  !   xswmin < xswinv < xswmix

  ! originally from "SCRATCH..."

  integer, dimension(:), allocatable :: iord,izonec
                           ! Helical flux ordering of zones, and
                           ! complementary zone with matching Helical flux
                           ! iord(j) -- zone bdy with j'th highest Helical
                           ! flux; izonec(j) -- zone # in which the other
                           ! itercept of helflx(iord(j),1) is located.

  real(fp), dimension(:,:), allocatable :: zxipre
                           ! min,max xi location of j'th Helical flux
                           ! inercept.

  integer, dimension(:), allocatable :: inzmix
                           ! inzmix(k) related to no. of pre-sawtooth zones
                           ! contributing material to k'th post-sawtooth
                           ! zone

  integer, dimension(:,:), allocatable :: izmix
  real(fp), dimension(:,:), allocatable :: zzmix
                           ! izmix(1:2,k),zzmix(1:2,k) give the pre-sawtooth
                           ! zone numbers and fractions specifying mixing into
                           ! post sawtooth zone k.

  integer :: imix          ! index to last zone involved in mixing
  integer :: imax          ! index of zone containing highest Helical flux
                           ! value (corresponds to q=1 surface).
  integer :: imin          ! index to zone containing inner bdy of q=1 island

  ! originally from "SCRATCH..." extension in sawnbi.for
  ! in the above arrays, one considers what pre-sawtooth zones map into
  ! a selected post-sawtooth zone; in these one considers what post-sawtooth
  ! zones are mapped to from a selected pre-sawtooth zone.

  integer, dimension(:), allocatable :: inzmxp
      !  INZMXP(J) = NUMBER OF POST SAWTOOTH ZONES MIXED INTO from zone J

  integer, dimension(:,:), allocatable :: izmixp
      !  IZMIXP(1:INZMXP(J),J) = THE ZONE NUMBERS OF THE MIXED-TO ZONES

  real(fp), dimension(:,:), allocatable :: zzmixp
      !  ZZMIXP(1:INZMXP(J),J) = the probability factor for each mixed-to zone.

end module kdsaw_module
