!-------------------------------------------------------------------
!  KDSAW_HDR
!  CALCULATE ORDER MAP FOR KADOMTSEV SAWTOOTH EVENT
!
!  MODIFICATIONS DMC NOV 1988 -- ADJUSTING DETAILS OF TREATMENT
!  OF HELICAL FLUX FUNCTION IN VICINITY OF ITS MAXIMUM (I.E. THE
!  Q=1 SURFACE)
!
SUBROUTINE KDSAW_HDR(ierr)
!
  use kdsaw_module                  ! used to be TRANSP COMMON
  implicit NONE
!
  integer, intent(out) :: ierr   ! error code, 0=OK
!
!  ierr=0: OK; ierr=99 -- error which indicates software bug!
!                         this is note expected to happen.
!
!---------------------------------------------------------------
!
!  INPUT-- TRANSP COMMON
!   HELFLX - THE PRE-SAWTOOTH HELICAL FLUX FUNCTION
!     WHICH HAS BEEN CHECKED IN THE SAWHCK SUBROUTINE
!          LOCAL COMMON /SCRATCH/
!   IMIX - THE LAST TRANSP ZONE INVOLVED IN SAWTOOTH MIXING
!   IMAX - THE TRANSP ZONE CONTAINING THE HIGHEST HELICAL FLUX VALUE
!
!   ORDER:
!    1<---REGION OF HIGHEST HELICAL FLUX (PRE-SAWTOOTH)
!         WHICH WILL BE MAPPED TO THE MAGNETIC AXIS (POST-SAWTOOTH)
!    2<---REGION OF NEXT HIGHEST HELICAL FLUX (PRE) WHICH WILL
!         BE MAPPED A LITTLE FURTHER FROM THE MAGNETIX AXIS (POST)
!   ...
!    N<---REGION OF HELICAL FLUX JUST .GT. ZERO (PRE) WHICH WILL BE
!         MAPPED TO NEAR THE OUTSIDE EDGE OF THE MIXING REGION (POST)
!
!  OUTPUTS  LOCAL COMMON /SCRATCH/  IORD AND IZONEC:
!    INDEX PTRS TO BIFURCATED REGIONS OF EQUAL HELICAL FLUX
!
!  REGIONS ARE DEFINED AS ZONES OF EQUAL PRE-SAWTOOTH HELICAL
!  FLUX.  THE (GREATER HELICAL FLUX) BDY OF THE JTH FLUX ZONE
!  CORRESPONDS TO TRANSP FLUX SURFACE IORD(J) WITH FLUX HELFLX(IORD(J)).
!  THESE REGIONS GENERALLY CONSIST OF TWO INTERVALS (EXCEPT FOR THE 1ST
!  REGION) AND THE COMPLIMENTARY INTERVAL HAS ITS BDY IN TRANSP ZONE
!  IZONEC(J).
!
!  THERE IS A SINGULARITY WHERE THE DERIVATIVE OF THE HELICAL FLUX
!  FCN GOES TO ZERO, I.E. AT THE MAXIMUM (HIGHEST ORDERED REGIONS).
!  THESE REGIONS MAP THE LARGEST AMOUNT OF TOROIDAL FLUX FOR THE
!  SMALLEST CHANGE IN HELICAL FLUX AND ARE THEREFORE GIVEN SPECIAL
!  TREATMENT (THE HELICAL FLUX FCN IS INTERPOLATED AS A PARABOLA
!  INSTEAD OF PIECE-WISE LINEAR)
!
!   ZXIPRE
!
!   ZXIPRE(1,J),ZXIPRE(2,J)
!  THIS ARRAY GIVES THE 2 XI (TOROIDAL FLUX LABEL) INTERCEPTS
!  CORRESPONDING TO THE SINGLE VALUE OF PRE-SAWTOOTH HELICAL FLUX
!  INDICATED BY IORD(J) (I.E. HELFLX(IORD(J)))
!
!-----------------------------------
!   THE OUTPUT OF THIS ROUTINE SUPPLIES THE INFORMATION NECESSARY
!  FOR THE SAWMIX ROUTINE TO CONSTRUCT THE POST SAWTOOTH HELICAL FLUX
!  AND IOTA(BAR) PROFILE AND THE PTCL/ENERGY MIXING FACTORS.
!
!   NOTE THAT THE RESULTING ORDERING OF TOROIDAL FLUX PAIRS VS.
!  HELICAL FLUX IS NOT IN EVEN STEPS OF HELICAL FLUX.  THE BDY'S ARE
!  DETERMINED BY WHERE THE TRANSP ZONE BDY'S HIT THE PRE-SAWTOOTH
!  HELICAL FLUX FUNCTION
!   ** IN FACT IF TWO TRANSP BDY'S HAPPEN TO HAVE EXACTLY THE SAME
!  PRE SAWTOOTH HELICAL FLUX, THIS ROUTINE RECORDS 2 CONSECUTIVE
!  IDENTICAL REGIONS (HELICAL FLUX STEPSIZE=0) *** CAUTION ***
!
!  LOCAL TO THIS ROUTINE:
!
  INTEGER IB(2),IC(2)
  integer jsign,imxp1,imxp2,inum,inp,istop,ip,iz,izp1,imxm1
  real(fp) :: zh0,zx0,zs,zdxi,zh1,zx1,zx2,zx1c,zh,zhp
  real(fp) :: za,zxmax,zx0c
!
!  EXECUTABLE CODE
!
!  1.  DEAL WITH THE REGIONS OF HIGHEST HELICAL FLUX
!      A PARABOLIC INTERPOLATION SCHEME IS USED.  THUS THE MAXIMUM
!      HELICAL FLUX IS NOT NECESSARILY DEFINED TO BE ON A TRANSP ZONE
!      BDY-- BUT IS TREATED AS THE PEAK OF THE PARABOLA DEFINED BY
!      3 PTS OF HELICAL FLUX AROUND THE MAXIMUM
!
!  SAWHCK CHECKS GUARANTEE THAT THE MAXIMUM HELICAL FLUX IS IN
!  THE INTERIOR BOUND AWAY FROM THE MAGNETIC AXIS
!
!  RECODED DMC NOV 1988 DUE TO DIFFICULTIES WITH FLUX PROFILES WITH
!  CERTAIN PROPERTIES NEAR THE Q=1 SURFACE
!
  ierr=0
!
  IF(HELFLX(IMAX+1).GT.HELFLX(IMAX)) THEN
!  (IF IMAX=LCENTR WE MUST COME HERE).
     JSIGN=-1
     IMXP1=IMAX+1
     IMXP2=IMAX+2
     !  (ACTUALLY - THE SLOPE HERE)
     ZH0=HELFLX(IMXP1)
     ZX0=XI(IMXP1,1)
     ZDXI=(XI(IMXP2,1)-XI(IMXP1,1))
     ZS=(ZH0-HELFLX(IMXP2))/ZDXI
     ZH1=HELFLX(IMAX)
     ZX1C=XI(IMAX,1)
  ELSE
     JSIGN=1
     IMXM1=IMAX-1
     IMXP1=IMAX+1
     ZH0=HELFLX(IMAX)
     ZX0=XI(IMAX,1)
     ZDXI=(XI(IMAX,1)-XI(IMXM1,1))
     ZS=(ZH0-HELFLX(IMXM1))/ZDXI
     ZH1=HELFLX(IMXP1)
     ZX1C=XI(IMXP1,1)
  ENDIF
!
  ZA=((ZH1-ZH0)-ZS*ZDXI)/(ZDXI*ZDXI)
!
!  PARABOLIC FUNCTION OF X FROM 0 TO ZDXI:
!   F(X)=ZA*X**2 + ZS*X + ZH0   GIVING F(ZDXI)=ZH1
!
!  F DEFINES INTERPOLATION OF HELICAL FLUX IN ZONE WHERE MAXIMUM (I.E.
!  Q=1) OCCURS.
!
!  LOCATE MAXIMUM
!
  ZXMAX=-ZS/(two*ZA)
  IF((ZXMAX.GT.half*ZDXI).AND.(ZXMAX.LT.0.500001_fp*ZDXI))        &   ! ???
         ZXMAX=half*ZDXI
!
  IF((ZXMAX.LE.0).OR.(ZXMAX.GT.half*ZDXI)) THEN
     WRITE(abs(NONLIN),9901)
 
 9901   FORMAT(' ?KDSAW_HDR -- ERROR IN TREATMENT OF HELICAL FLUX FUNCTION')
     ierr=99
     return
  ENDIF
  !
  !  FIRST ZONE
  IORD(1)=IMAX-(JSIGN-1)/two
  IZONEC(1)=IMAX
 
  ZX0C=ZX0+JSIGN*two*ZXMAX
  IF(JSIGN.EQ.1) THEN
     ZXIPRE(1,1)=ZX0
     ZXIPRE(2,1)=ZX0C
  ELSE
     ZXIPRE(1,1)=ZX0C
     ZXIPRE(2,1)=ZX0
  ENDIF
  !
  !-----------------------------------
  !  DMC *BUGFIX*  XSWINV = INVERSION LOCN MUST BE PLACED IN MIDDLE OF
  !   THIS FIRST RECONECTION ZONE
  !
  XSWINV=(ZXIPRE(1,1)+ZXIPRE(2,1))/two
  !
  !  2.  USING LINEAR INTERPOLATION OF THE HELICAL FLUX FUNCTION, DEAL
  !      WITH ALL REMAINING REGIONS
  !
  INUM=1
  ISTOP=IMIX-LCENTR+1
  !
 10   CONTINUE
  !  LOOP, EXIT CONDITION
  IF(INUM.EQ.ISTOP) GO TO 100
  !
  INP=INUM
  ZHP=HELFLX(IORD(INP))
 
  INUM=INUM+1
 
  IF(IZONEC(INP).LT.IORD(INP)) THEN
     !  FIND PTRS TO EACH OF THE TWO POSSIBLE NEXT BDY'S
     !  SCANNING DEPENDS ON "WHICH SIDE" OF THE HELICAL FLUX .GT. 0
     !  REGION THE LAST BDY WAS ON.
     !
     !  2 SIDES ARE SIDE 1:  D/D(XI)[HELFLX] .GT. 0 NEAR THE MAGNETIC AXIS
     !              SIDE 2:  D/D(XI)[HELFLX] .LT. 0 AWAY FROM MAGNETIC AXIS
     !  PREVIOUS BDY ON SIDE 2
     IB(1)=IZONEC(INP)
     IB(2)=IORD(INP)+1
     !  ALSO DEFINE OPPOSITE SIDE INTERSECTION PTRS
     IC(1)=IORD(INP)
     IC(2)=IZONEC(INP)
  ELSE
     !  PREVIOUS BDY ON SIDE 1
     IB(1)=IORD(INP)-1
     IB(2)=IZONEC(INP)+1
     IC(1)=IZONEC(INP)
     IC(2)=IORD(INP)-1
  ENDIF
  !  THE NEXT ORDER BDY IS THE ONE WITH THE GREATER HELICAL FLUX
  IF(HELFLX(IB(1)).GT.HELFLX(IB(2))) THEN
     IP=1
  ELSE IF(HELFLX(IB(1)).LT.HELFLX(IB(2))) THEN
     IP=2
   ELSE
     !  REMOTE CHANCE: HELICAL FLUX EQUALITY OF BDY'S:  SPECIAL HANDLING
     IF(HELFLX(IB(1)).EQ.ZHP) THEN
        IP=2
     ELSE
        IP=1
     ENDIF
  ENDIF
 
  IORD(INUM)=IB(IP)
  IZONEC(INUM)=IC(IP)
  ZH=HELFLX(IORD(INUM))
  ZX1=XI(IORD(INUM),1)
  IZ=IZONEC(INUM)
  IF(IZ.NE.IMAX) THEN
     !  STRAIGHT LINEAR INTERPOLATION
     IZP1=IZ+1
     ZX2=XI(IZ,1)+(XI(IZP1,1)-XI(IZ,1))*                     &
         (ZH-HELFLX(IZ))/(HELFLX(IZP1)-HELFLX(IZ))
  ELSE
     !  INTERCEPT IN ZONE SHARED WITH HELICAL FLUX MAXIMUM.
     ZX2=ZX0C+(ZX1C-ZX0C)*(ZH-ZH0)/(ZH1-ZH0)
  ENDIF
!
  ZXIPRE(1,INUM)=MIN(ZX1,ZX2)
  ZXIPRE(2,INUM)=MAX(ZX1,ZX2)
!
  GO TO 10
!
!  EXIT
!
 100   CONTINUE
!
  return    
end subroutine kdsaw_hdr
