subroutine kdsaw_find_r4(zselect,izone_in,izone_out,zfrac)

  ! REAL interface to kdsaw_find

  implicit NONE

  real,intent(in) :: zselect  ! a number btw 0 and 1
  integer,intent(in) :: izone_in  ! pre-sawtooth originating zone
  integer, intent(out) :: izone_out  ! post-sawtooth zone selected
  real, intent(out) :: zfrac  ! fractional location within zone

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

  real*8 zselect_r8, zfrac_r8

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

  zselect_r8 = zselect

  call kdsaw_find(zselect_r8,izone_in,izone_out,zfrac_r8)

  zfrac = zfrac_r8

end subroutine kdsaw_find_r4

subroutine kdsaw_find(zselect,izone_in,izone_out,zfrac)

  use kdsaw_module
  implicit NONE

  ! given an input zone and a (random) number between 0 and 1, select
  ! an output zone, according to the Kadomtsev redistribution mapping

  real(fp),intent(in) :: zselect  ! a number btw 0 and 1
  integer,intent(in) :: izone_in  ! pre-sawtooth originating zone
  integer, intent(out) :: izone_out  ! post-sawtooth destination zone selected
  real(fp), intent(out) :: zfrac  ! fractional location within zone

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

  real(fp) :: zsum,zsump,zsel_cor
  integer i,izi,izo

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

  izi=izone_in+lcentr-1

  if(inzmxp(izi).eq.0) then
     write(abs(nonlin),*) ' %kdsaw_find: warning: izone_in = ',izone_in, &
          ' not in mixing region.'
     izone_out = izone_in
     zfrac=half
  else if((zselect.lt.0).or.(zselect.gt.1)) then
     write(abs(nonlin),*) ' %kdsaw_find: warning: zselect = ',zselect, &
          ' not between 0 and 1.'
     izone_out = izone_in
     zfrac=half 
  else
     !  search adjusted due to fact that last zone is only partially
     !  included in mixing region

     if(izi.lt.imix) then
        zsel_cor = zselect
     else
        zsel_cor=zero
        do i=1,inzmxp(izi)
           izo=izmixp(i,izi)
           if(izo.lt.imix) then
              zsel_cor=zsel_cor+zzmixp(i,izi)
           else
              zsel_cor=zsel_cor+zzmixp(i,izi)*zfrmix
           endif
        enddo
        zsel_cor = zsel_cor*zselect
     endif

     zsum=zero
     do i=1,inzmxp(izi)
        izo=izmixp(i,izi)
        zsump=zsum
        if(min(izi,izo).lt.imix) then
           zsum=zsum+zzmixp(i,izi)
        else
           zsum=zsum+zzmixp(i,izi)*zfrmix
        endif
        if((zsump.le.zsel_cor).and.(zsel_cor.le.zsum)) exit
     enddo

     zfrac=(zsel_cor-zsump)/(zsum-zsump)
     if(izo.eq.imix) then
        zfrac=zfrac*(xswmix-xi(imix,1))/(xi(imix+1,1)-xi(imix,1))
     endif

     izone_out=izo-lcentr+1

  endif

end subroutine kdsaw_find
