!-----------------------------------------------------------------
!  KDSAW_HCK
!  CHECK THE PRE-SAWTOOTH HELICAL FLUX FUNCTION.
!    REQUIRE:
!  1.  THE DERIVATIVE OF THE HELICAL FLUX FCN BE .GT. 0.0 ON AXIS
!      (EQUIVALENT TO Q<1 ON AXIS, EXPECTED FOR A SAWTOOTH)
!  2.  EXISTANCE OF ONE AND ONLY ONE POINT WHERE THE DERIVATIVE OF
!      THE HELICAL FLUX FCN IS ZERO (EQUIVALENT TO EXISTANCE AND
!      UNIQUENESS OF A Q=1 SURFACE)
!  3.  THE HELICAL FLUX FCN BE .LT. 0.0 AT THE PLASMA BDY.  I.E.
!      THE KADOMTSEV MIXING REGION MUST NOT ENCOMPASS THE ENTIRE
!      PLASMA.  NOTE THAT BY TRANSP CONVENTION HELICAL FLUX =0
!      AT THE MAGNETIC AXIS.
!
!  ASSUMING THESE REQUIREMENTS ARE MET,
!
!    EVALUATE:
!      IMIX   INDEX TO THE LAST TRANSP ZONE AND ZONE BDY INSIDE
!             THE KADOMTSEV MIXING REGION
!      IMAX   INDEX TO THE ZONE IN WHICH THE Q=1 SURFACE LIES
!             (ZONE OF MAXIMUM HELICAL FLUX)
!      IMIN   INDEX to zone containing inner bdy of q=1 island; 
!             for Kadomtsev this is LCENTR always; for Porcelli it
!             can be further out.
!
!      XSWMIX RELATIVE FLUX LABEL (XI) CORRESPONDING TO BDY OF
!             KADOMTSEV MIXING REGION
!      XSWINV RELATIVE FLUX LABEL (XI) CORRESPONDING TO Q=1 SURFACE
!             I.E. SAWTOOTH INVERSION LAYER
!      XSWMIN RELATIVE FLUX LABEL (XI) CORRESPONDING TO inner bdy
!             of q=1 island.
!
!      ZVOLAX volume of axial region enclosed by q=1 island region
!      ZVOLQ1 volume of q=1 island region
!      ZVOLMIX  volume of total mixing region = ZVOLAX+ZVOLQ1
!
!      ZFRMIX fraction of volume of last zone that is in mixing region
!      ZFRMIN fraction of volume of last axial region zone that is in
!             axial region
!
!------------------------------------------------------------
!
SUBROUTINE KDSAW_HCK(IERR)
!
  use kdsaw_module                  ! used to be TRANSP COMMON
  implicit NONE
!
  integer,intent(out) :: ierr  ! set only if q<1 all the way to outer bdy
!
!------------------------
  integer j,jp1,ionlin
  real(fp) :: zib1,zib2,zx1,zx2,zhf1,zhf2,zdv
!-------------------------------------------------------------------
!  EXECUTABLE CODE
!
  IERR=0
!
  ionlin=abs(nonlin)
!
  call clear                ! initially assume no mixing...
!
!  REQUIREMENT 1. TANTAMOUNT TO 1ST HELICAL FLUX VALUE .GT. 0.0
!
  IF(HELFLX(LCP1).LE.zero) THEN
     RETURN                         ! no mixing: q>1 on axis
  ENDIF
!
!  REQUIREMENT 3. THE HELICAL FLUX AT THE EDGE MUST BE NEGATIVE
!
  IF(HELFLX(LEP1).GE.zero) THEN
     write(ionlin,*) ' ?kdsaw_hck:  q=1 too close to boundary!'
     IERR=1
     RETURN
  ENDIF
!
!  REQUIREMENT 1+3 SATISFIED IMPLIES EXISTANCE OF Q=1 SURFACE.
!  NOW CHECK UNIQUENESS WHILE EVALUATING ROUTINE OUTPUTS.
!  SCAN THE IOTA(BAR) PROFILE
!   NOTE RELATION D(HELICAL FLUX)/D(XI)
!         = ((IOTA(BAR)-1)*XI*TFLUX)/PI
!   AND IOTA(BAR)=1/Q
!
  DO J=LCENTR,LEDGE
     JP1=J+1
     ZIB1=XIOTBP(J,1)
     ZIB2=XIOTBP(JP1,1)
     ZX1=XI(J,1)
     ZX2=XI(JP1,1)
     ZHF1=HELFLX(J)
     ZHF2=HELFLX(JP1)
     !  UNIQUENESS CHECK:  ONCE .GE. 1 Q MUST NEVER GET .LE. 1 AGAIN
     IF(IMAX.GT.zero) THEN
        IF(ZIB2.GE.one) then
           IERR=2
           write(ionlin,*) ' ?kdsaw_hck:  q=1 surface *not unique*!'
           write(ionlin,*) '  ...sawtooth mixing suppressed.'
           call clear
           return
        endif
     endif
     !  LOOK FOR OCCURANCE OF Q=1 SURFACE
     IF((ZIB1.GT.one).AND.(ZIB2.LE.one)) THEN
        !  Q=1 SURFACE FOUND
        IMAX=J
        XSWINV=ZX1+(ZX2-ZX1)*(one-ZIB1)/(ZIB2-ZIB1)
     ENDIF
     !  LOOK FOR END OF MIXING REGION = 2ND SURFACE WHERE HELEICAL FLUX = 0
     if(mix_method.le.2) then
        IF((ZHF1.GT.zero).AND.(ZHF2.LE.zero)) THEN
           ! KADOMTSEV Mixing region -- Helical flux crosses zero at outer edge
           ! END OF MIXING REGION FOUND
           IMIX=J
           XSWMIX=ZX1-(ZX2-ZX1)*ZHF1/(ZHF2-ZHF1)
           zfrmix=(XSWMIX**two-xi(j,1)**two)/(xi(jp1,1)**two-xi(j,1)**two)
           ZVOLQ1=ZVOLQ1+DVOL(J,2)*zfrmix
        ELSE if(zhf2.gt.zero) then
           ZVOLQ1=ZVOLQ1+DVOL(J,2)
        ENDIF
     endif
     !  LOOP END
  enddo
 
  if(mix_method.le.2) then
     !  Kadomtsev: no axial region, q=1 island extends all the way to axis.
     zfrmin=zero
     XSWMIN=zero
     ZVOLAX=zero
     ZVOLMIX=ZVOLQ1
     imin=lcentr
  else
     !  Porcelli method
     XSWMIN=f_porcelli*XSWINV
     do j=lcentr,ledge
        jp1=j+1
        if(xi(jp1,1).lt.XSWMIN) then
           !  add to volume of axial region
           zvolax=zvolax+dvol(j,2)
           cycle
        else if(imin.eq.0) then
           !  inside bdy of q=1 island found...
           !  end of axial region, start of q=1 island region
           imin=j
           if(imin.eq.imax) then
              ! special handling: q=1 island is single zone...
              imix=imin
              zfrmin=zero
              zfrmix=one
              zvolq1=dvol(j,2)
              XSWMIN=xi(j,1)
              XSWMIX=xi(jp1,1)
              exit
           else
              zfrmin=(XSWMIN**two-xi(j,1)**two)/(xi(jp1,1)**two-xi(j,1)**two)
              zdv=dvol(j,2)*zfrmin
              zvolax=zvolax+zdv
              zvolq1=dvol(j,2)-zdv
              zhf1=helflx(j)+(XSWMIN-xi(j,1))*                  &
                   (helflx(jp1)-helflx(j))/(xi(jp1,1)-xi(j,1))
           endif
        else
           !  looking for outer q=1 island boundary by matching
           !  helical flux at inner bdy (Porcelli method)
           if((helflx(j).ge.zhf1).and.(helflx(jp1).le.zhf1)) then
              imix=j
              if(helflx(j).eq.helflx(jp1)) then
                 xswmix=xi(j,2)
              else
                 xswmix=xi(j,1)+(xi(j,2)-xi(j,1))*           &   
                      (zhf1-helflx(j))/(helflx(jp1)-helflx(j))
              endif
              xswmix=max(xswinv,xswmix)
              zfrmix=(XSWMIX**two-xi(j,1)**two)/                &
                     (xi(jp1,1)**two-xi(j,1)**two)
              zdv = dvol(j,2)*zfrmix
              zvolq1=zvolq1+zdv
              exit
           else
              zvolq1=zvolq1+dvol(j,2)
           endif
        endif
     enddo
     zvolmix=zvolq1+zvolax
  endif
!
!  EXIT
!
  RETURN
!
!----------------------------
!
contains
!
subroutine clear
!
  imix=0
  imax=0
  imin=0
!
  xswmix=zero
  xswinv=zero
  xswmin=zero
!
  zvolq1=zero
  zvolax=zero
  zvolmix=zero
!
  zfrmix=zero
  zfrmin=zero
!
end subroutine clear

END subroutine kdsaw_hck
