C--------------------------------------------------------------------
C
      SUBROUTINE KDSAW_MIX(ierr)
C
      use kdsaw_module                  ! used to be TRANSP COMMON
      implicit NONE
C     
      integer, intent(out) :: ierr      ! error code, 0=OK
C
C  ierr=0: OK; ierr=99 -- error which indicates software bug!
C                         this is not expected to happen.
C
C---------------------------------------------------------------
C
C  KDSAW_MIX - DMC 27 JULY 1984
C   USING THE INFORMATION PREPARED IN SUBROUTINE SAWHDR,
C  CALCULATE THE KADOMTSEV SAWTOOTH MIXING FACTORS.
C
C   BASED ON EQUALITY OF HELICAL FLUX, BIFURCATED REGIONS OF
C  TOROIDAL FLUX WERE MAPPED OUT IN THE SAWHDR ROUTINE.  IN THIS
C  ROUTINE THESE REGIONS ARE COMBINED AND RE-ORDERED TO GET THE
C  POST SAWTOOTH IOTA(BAR) (EQUIV. PLASMA CURRENT) AND THERMAL
C  PRESSURE PROFILES.  THE RE-ORDERING CONSERVES TOROIDAL FLUX,
C  I.E. THE FLUX IN POST-SAWTOOTH ZONE J EQUALS THE SUM OF
C  FLUXES IN PRE-SAWTOOTH ZONES K1, K2, ... KN FROM WHENCE THE
C  CURRENT AND MATERIAL IN ZONE J CAME.  (TOROIDAL FLUX CONSERVATION
C  EQUIVALENT TO PRESERVING XI**2 = INTEGRAL (XI*DXI) WHERE XI IS
C  TRANSP NORMALIZED TOROIDAL FLUX COORDINATE)
C
C-------------------------------------------------------------------
C  DMC November 2005 -- mix_method control added.
C  If mix_method.eq.1 the old Kadomtsev helical flux matching method
C    is in effect.
C  Otherwise:
C    mix_method=2 or mix_method=4: even mixing throughout region
C    mix_method=3: even mix in axial region AND separate even mix in q=1
C                  island region.
C
C-------------------------------------------------------------------
C  SUBROUTINE USES SAWHDR OUTPUTS AS INPUTS + THIS INFO FROM
C   TRANSP COMMON:
C   XI  THE RELATIVE TOROIDAL FLUX AT ZONE BDYS
C   DVOL  ZONE VOLUMES-- USED IN FIGURING THE WEIGHTED MIXING FACTORS
C   XIOTBP  THE PRE-SAWTOOTH IOTA(BAR) PROFILE
C
C  OUTPUT-- TO LOCAL COMMON BLOCK:
C   FOR POST-SAWTOOTH ZONE J,
C
C   INZMIX(J)  NUMBER OF PIECES OF PRE-SAWTOOTH ZONES MIXED IN
C              TO FORM THIS POST-SAWTOOTH ZONE
C   IZMIX(....,J)  THE ZONE NUMBERS OF THE MIXED IN ZONES
C   ZZMIX(....,J)  THE MIXING WEIGHTING FACTORS FOR ZONE J
C           TO TRANSP COMMON:
C   XIOTB   AT ZONE BDYS-- THE POST SAWTOOTH IOTA(BAR) PROFILE
C
C  **MODIFIED DMC OCT 1986** USE NAMELIST CONTROL XSWID2
C
C  IF XSWID2=0 (DEFAULT) STANDARD MODEL
C  IF XSWID2=1 *DO NOT CHANGE THE IOTA(BAR) PROFILE*
C  IF XSWID2 BTW 0 AND 1 DO WEIGHTED AVERAGE OF THESE EXTREMA
C
C-------------------------------------------------------------------
C  STATEMENT FUNCTIONS--
C
      integer k
      real(fp) :: zib,zibps,zib0,zib1,zib2,zx,zdx,zsib,zw
C
C  INTERPOLATE IOTA(BAR) PROFILE
C
      ZIB(K,ZX)=XIOTBP(K,1)+(XIOTBP(K+1,1)-XIOTBP(K,1))*
     >            (ZX-XI(K,1))/(XI(K+1,1)-XI(K,1))
C
C  CALCULATE POST-SAWTOOTH IOTA(BAR) FROM TWO PRE-SAWTOOTH IOTA(BAR)
C  VALUES AT THE SAME HELICAL FLUX
C
      ZIBPS(ZIB1,ZIB2)= one - (ZIB1-one)*(ZIB2-one)/(ZIB2-ZIB1)
C
C  local variables
C
      integer ip,ipstop,ix1,ix2,iz
      real(fp) :: zx1,zx2,zx1p,zx2p,zf,zfp,zfps,zfpsp
      integer j,ii,ij,jj,jp1
      real(fp) :: zintrp,zx1i,zx2i
C
      real(fp) :: ziotb(lep1)
C
C------------------------------------------------------------
C
      ierr=0
C
      if(mix_method.gt.2) go to 500
C
C  EXECUTABLE CODE
C  INITIATE IOTA(BAR)=1 ON AXIS
C
      XIOTB(LCENTR,1)=1
C
C  INITIAL PRE-SAWTOOTH BDY:  THE INVERSION LAYER
C
      IP=0
      IPSTOP=IMIX-LCENTR+1
C
      IX1=IMAX
      IX2=IMAX
C
      ZX1=XSWINV
      ZX1P=ZX1
      ZX2=XSWINV
      ZX2P=ZX2
      ZF=zero
      ZFP=ZF
C
      ZFPS=zero
      ZFPSP=ZFPS
C
C  LOOP THRU (POST-SAWTOOTH) ZONES
C
      DO 100 J=LCENTR,IMIX
C
        JJ=J
        JP1=J+1
C  PREVIOUS ZONE, AND THIS ZONE'S
C  ENCLOSED (RELATIVE) TOROIDAL FLUX POST SAWTOOTH
        ZFPSP=ZFPS
        ZFPS=XI(JP1,1)**2
C
C  CLEAR MIX COUNTER
        INZMIX(J)=0
C  GET ANY LEFT OVERS FROM THE LAST PRE-SAWTOOTH ZONE PAIR
CX	  write(nonlin,'('' zf,zfpsp = '',2(1x,1pe11.4))') zf,zfpsp
        CALL KDSAWMF2(JJ,IX1,IX2,ZX1,ZX1P,ZX2,ZX2P,ZF,ZFP,ZFPS,ZFPSP)
CX	  write(nonlin,'('' 1)INZMIX('',i2,'')='',i2)') jj,inzmix(jj)
C
C
C  LOOK FOR PRESAWTOOTH REGIONS WHICH ENCLOSE THIS MUCH FLUX
C   OR END OF SAWTOOTH MIXING REGION
C
        DO WHILE ((ZF.LT.ZFPS).AND.(IP.NE.IPSTOP))
C
          ZX1P=ZX1
          ZX2P=ZX2
          ZFP=ZF
C
C  FIND NEXT PRE-SAWTOOTH BDY PAIR
          DO WHILE (ZF.LE.ZFP)
            IP=IP+1
            ZX1=ZXIPRE(1,IP)
            ZX2=ZXIPRE(2,IP)
            ZF=(ZX2-ZX1)*(ZX2+ZX1)
C  BEWARE OF POSSIBLE BDY PAIR DUPLICATION (OCCURS IF TWO PRE-SAWTOOTH
C  ZONE BDYS CORRESPONDED TO EXACTLY EQUAL TOROIDAL FLUX
          END DO
C
C  ZONE INDEX PTRS
          IF(IORD(IP).LE.IZONEC(IP)) THEN
            IX1=IORD(IP)
            IX2=IZONEC(IP)
          ELSE
            IX1=IZONEC(IP)
            IX2=IORD(IP)-1
          ENDIF
C  ADD MIXIN PIECES FROM THESE ZONES
          CALL KDSAWMF2(JJ,IX1,IX2,ZX1,ZX1P,ZX2,ZX2P,ZF,ZFP,ZFPS,ZFPSP)
CX	    write(nonlin,'('' 2)INZMIX('',i2,'')='',i2)') jj,inzmix(jj)
C  CHECK FOR ANOTHER PAIR OF PIECES
        END DO
C
C  DEAL WITH THE LAST ZONE-- INCLUDE IN THE UNMIXED PORTION
C
        IF((IP.EQ.IPSTOP).AND.(ZF.LT.ZFPS)) THEN
          IZ=INZMIX(J)+1
          INZMIX(J)=IZ
          IZMIX(IZ,J)=IMIX
          ZZMIX(IZ,J)=(ZFPS-ZF)/(ZFPS-ZFPSP)
        ELSE IF(ZF.GE.ZFPS) THEN
C
C  CALCULATE THE NEW BDY IOTA(BAR)
C
          ZINTRP=(ZFPS-ZFP)/(ZF-ZFP)
          ZX1I=ZX1P+ZINTRP*(ZX1-ZX1P)
          ZX2I=ZX2P+ZINTRP*(ZX2-ZX2P)
          ZIB1=ZIB(IX1,ZX1I)
          ZIB2=ZIB(IX2,ZX2I)
          XIOTB(JP1,1)=ZIBPS(ZIB1,ZIB2)
C
        ENDIF
C
C  BUG CHECK
C
        IF(INZMIX(J).LT.2) THEN
          WRITE(NONLIN,9901) J
 9901     FORMAT(' ?KDSAW_MIX BUGCHECK:  INZMIX(',I2,').LT.2'/)
          DO 99 II=1,J
            WRITE(NONLIN,9902) II,INZMIX(II)
 9902       FORMAT(' ...INZMIX(',I2,') = ',I2)
            DO 98 IJ=1,INZMIX(II)
              WRITE(NONLIN,9903) IJ,II,IZMIX(IJ,II),IJ,II,ZZMIX(IJ,II)
 9903         FORMAT(' .....IZMIX(',I2,',',I2,')=',I2/
     >                 ' .....ZZMIX(',I2,',',I2,')=',1PE11.4)
 98         CONTINUE
 99       CONTINUE
          ierr=99
        ENDIF
C
 100  CONTINUE
C
C------------------------------------------------------------
C  dmc Nov 2005 -- implement new mixing options
C
 500  continue
      if((ierr.eq.0).and.(mix_method.gt.1)) then

         call kdsaw_mix2(ierr)  ! see below...

      endif
C
C------------------------------------------------------------
C  option to smooth the q profile transitions
C
      if(dxi_qjump.gt.ZERO) then
         if(xswmin.gt.ZERO) then
C
C  smooth at inner bdy 
C
            zib0=xiotb(lcentr,1)
            zib1=ONE
            do j=lcentr,lep1
               if(xi(j,1).ge.xswmix) exit ! done if reach mixing bdy
               zx = xi(j,1)-xswmin
               if(zx.ge.dxi_qjump) exit ! done if beyond xswmin+dxi_qjump
               if(zx.lt.ZERO) then
                  if(zx.le.-dxi_qjump) cycle ! out of smoothing range low
                  zw = HALF + zx/dxi_qjump + 
     >                 zx*zx/(2*dxi_qjump*dxi_qjump)
               else
                  zw = HALF + zx/dxi_qjump -
     >                 zx*zx/(2*dxi_qjump*dxi_qjump)
               endif
               xiotb(j,1)=(ONE-zw)*zib0 + zw*zib1
            enddo
         endif
C
C  smooth at outer bdy
C
         zib2=xiotbp(imix+1,1)
         zx2=xi(imix+1,1)
         zsib=(zib2-xiotbp(imix,1))/(zx2-xi(imix,1))

         ziotb(lcentr:lep1)=xiotb(lcentr:lep1,1)

         do j=lcentr,lep1
            zx = xi(j,1) - xswmix
            if(zx.ge.dxi_qjump) exit ! done if beyond xswmix+dxi_qjump
            if(zx.lt.ZERO) then
               if(zx.le.-dxi_qjump) cycle ! out of smoothing range low
               zw = HALF + zx/dxi_qjump + 
     >              zx*zx/(2*dxi_qjump*dxi_qjump)
            else
               zw = HALF + zx/dxi_qjump -
     >              zx*zx/(2*dxi_qjump*dxi_qjump)
            endif
            if(j.le.imix) then
               zib0=xiotb(j,1)
            else if(imix.eq.lcentr) then
               zib0=xiotb(imix,1)
            else
               zib0=xiotb(imix,1)+(xi(j,1)-xi(imix,1))*
     >              (xiotb(imix,1)-xiotb(imix-1,1))/
     >              (xi(imix,1)-xi(imix-1,1))
            endif
            if(j.gt.imix) then
               zib1=xiotb(j,1)
            else
               zdx=zx2-xi(j,1)
               zib1=zib2-zsib*zdx
            endif
            ziotb(j)=(ONE-zw)*zib0 + zw*zib1
         enddo

         xiotb(lcentr:lep1,1)=ziotb(lcentr:lep1)
      endif
C    
      RETURN
C
      END
C------------------------------------------------------------
C
      subroutine kdsaw_mix2(ierr)
C
      use kdsaw_module                  ! used to be TRANSP COMMON
      implicit NONE
C     
      integer, intent(out) :: ierr      ! error code, 0=OK
C
C  ierr=0: OK; ierr=99 -- error which indicates software bug!
C                         this is not expected to happen.
C
C  mix_method=2 -- iota(bar) profile already computed
C  mix_method=3 or 4 -- iota(bar) profile still needs to be computed.
C
C  mix_method=2 or 4 --
C          compute "full random redistribution over entire mixing region"
C  mix_method=3 -- 
C          compute random redistribution separately within each of the
C          two sub-regions: the q=1 island region and the enclosed axial
C          region.
C
      real*8 :: zib0,zfrac,zfracmx
      integer :: j,jj,ict
C
C-------------------------------------
C
      ierr=0
      if(mix_method.eq.1) then
         write(nonlin,*) ' ?? kdsaw_mix2 -- invalid call, mix_method=1'
         ierr=99
         return
      endif
C
C  compute iota(bar) profile if necessary
C
      if(mix_method.gt.2) then

         zib0=xiotbp(imin,1)+(xswmin-xi(imin,1))*
     >        (xiotbp(imin+1,1)-xiotbp(imin,1))/
     >        (xi(imin+1,1)-xi(imin,1))

         do j=lcentr,lep1
            if(xi(j,1).le.xswmin) then
               xiotb(j,1)=zib0
            else if(xi(j,1).le.xswmix) then
               xiotb(j,1)=ONE
            else
               xiotb(j,1)=xiotbp(j,1)
            endif
         enddo
      endif
C
      if((mix_method.eq.2).or.(mix_method.eq.4)) then
C
C  compute even mixing over entire sawtooth region
C
         do j=lcentr,imix
            if(j.lt.imix) then
               zfrac=dvol(j,2)/zvolmix
            else
               zfrac=zfrmix*dvol(j,2)/zvolmix
            endif
            inzmix(j)=imix-lcentr+1
            ict=0
            do jj=lcentr,imix-1
               ict=ict+1
               izmix(ict,j)=jj
               zzmix(ict,j)=zfrac
            enddo
            ict=ict+1
            izmix(ict,j)=imix
            zzmix(ict,j)=zfrmix*zfrac  ! fraction zfrmix of zone is sent
            !  and the part of the last zone not in the mixing:
            if(j.eq.imix) zzmix(ict,j)=zzmix(ict,j)+(1-zfrmix)
         enddo
C
      else
C
C  compute even mixing within axial region, separated from even mixing
C    within q=1 island region
C
C  special case: q=1 island restricted to a single zone:
C    cf kdsaw_hck, in this case, the q=1 island is treated as coincident
C    with the zone contain the q=1 surface; the axial region is all zones
C    inside of but not including this zone.
C
         if(imin.eq.imix) then
            inzmix(imin)=1
            izmix(1,imin)=imin
            zzmix(1,imin)=1
            do j=lcentr,imin-1
               zfrac=dvol(j,2)/zvolax
               inzmix(j)=imin-lcentr
               ict=0
               do jj=lcentr,imin-1
                  ict=ict+1
                  izmix(ict,j)=jj
                  zzmix(ict,j)=zfrac
               enddo
            enddo
C
         else
C
C  normal case: q=1 inner and outer boundaries in different zones...
C
C  compute even mixing within axial region
C
C  this loop:  all zones entirely within axial region
C
            do j=lcentr,imin-1
               zfrac=dvol(j,2)/zvolax
               inzmix(j)=imin-lcentr+1
               ict=0
               do jj=lcentr,imin-1
                  ict=ict+1
                  izmix(ict,j)=jj
                  zzmix(ict,j)=zfrac
               enddo
               ict=ict+1
               izmix(ict,j)=imin
               zzmix(ict,j)=zfrmin*zfrac
            enddo
C
C  here: the zone (imin) split btw axial region and q=1 island region
C
            j=imin
            zfrac=zfrmin*dvol(j,2)/zvolax
            inzmix(j)=imix-lcentr+1 ! include contributions from all zones
            ict=0
            do jj=lcentr,imin-1
               ict=ict+1
               izmix(ict,j)=jj
               zzmix(ict,j)=zfrac
            enddo
            ict=ict+1
            izmix(ict,j)=imin
            zzmix(ict,j)=zfrmin*zfrac
C
C  axial region to (imin) done; add in q=1 region contribution
            zfrac=(1-zfrmin)*dvol(j,2)/zvolq1
            zzmix(ict,j)=zzmix(ict,j)+(1-zfrmin)*zfrac
C
C  add in rest of q=1 region contributions to (imin) boundary zone
C
            do jj=imin+1,imix-1
               ict=ict+1
               izmix(ict,j)=jj
               zzmix(ict,j)=zfrac
            enddo
            ict=ict+1
            izmix(ict,j)=imix
            zzmix(ict,j)=zfrmix*zfrac
C
C  this loop computes mixing of zones entirely within q=1 island region
C
            do j=imin+1,imix
               if(j.lt.imix) then
                  zfrac=dvol(j,2)/zvolq1
               else
                  zfrac=zfrmix*dvol(j,2)/zvolq1
               endif
               inzmix(j)=imix-imin+1
               ict=1
               izmix(ict,j)=imin
               zzmix(ict,j)=zfrac*(1-zfrmin)
               do jj=imin+1,imix-1
                  ict=ict+1
                  izmix(ict,j)=jj
                  zzmix(ict,j)=zfrac
               enddo
               ict=ict+1
               izmix(ict,j)=imix
               zzmix(ict,j)=zfrmix*zfrac ! fraction zfrmix of zone is sent
               !  and the part of the last zone not in the mixing:
               if(j.eq.imix) zzmix(ict,j)=zzmix(ict,j)+(1-zfrmix)
            enddo
         endif                  ! single vs. multi-zone q=1 island
C
      endif                     ! mix all vs. mix axial, q=1 regions separately
C
      return
      end
