#include "fpreproc/f77_dcomplx.h"
      subroutine rfgbmmx(zthBmin,zBmin,zthBmax,zBmax)
c
c dmc at JET -- Aug 1995.  For FPPMOD with up-down asymmetry.
c   return theta(Bmin) and Bmin, theta(Bmax) and Bmax.
c
c   flux surface zxi and R*Bt in rfgcom COMMON have already been set ***
c
c *B in gauss*
c *internal use only*
C============
C idecl:  explicitize implicit INTEGER declarations:
!
      IMPLICIT NONE
      REAL*8, PARAMETER :: HALF = 0.5d0
      REAL*8, PARAMETER :: ONE = 1.0d0
      REAL*8, PARAMETER :: TWO = 2.0d0
!
      INTEGER isrch,imin,imax
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zthbmin,zbmin,zthbmax,zbmax,zpi,zpio2,ztwopi,zth0,
     > zth1,rfgbinv,zthmin0,zbmin0,zthmin1,zbmin1,zthmax0,zbmax0,
     > zthmax1,zbmax1
C============
      REAL*8 zBrequ ! requested B field value.
      REAL*8 zxi    ! on this flux surface (normalized rho/rho_lim)
      REAL*8 zRB    ! R*Bt for rfgbcom
      common /rfgbcom/ zxi,zBrequ, zRB

      parameter (ISRCH=12)
      REAL*8 ztha(ISRCH),zBa(ISRCH)
c
c  this source file includes subroutine rfgbmmx, rfgbmmx0, and rfgbmmx1.
c  subfunction rfgbinv is used to evaluate B(theta).
c
c--------------------
c
      zpi=abs(acos(-ONE))
      zpio2=HALF*zpi
      ztwopi=TWO*zpi
c
c  COMMON:
      zBrequ=0.0  ! so rfgbinv returns B

c  initial rough search

      zth0=-zpio2
      zth1=zth0+ztwopi

      ztha(1)=zth0
      zBa(1)=rfgbinv(zth0)
      ztha(ISRCH)=zth1
      zBa(ISRCH)=rfgbinv(zth1)
      call rfgbmmx0(ztha,zBa,ISRCH,imin,imax,0)

      zthmin0=ztha(imin-1)
      zBmin0=zBa(imin-1)
      zthmin1=ztha(imin+1)
      zBmin1=zBa(imin+1)

      zthmax0=ztha(imax-1)
      zBmax0=zBa(imax-1)
      zthmax1=ztha(imax+1)
      zBmax1=zBa(imax+1)

c  final searches for min and max B and thetas

      ztha(1)=zthmin0
      zBa(1)=zBmin0
      ztha(ISRCH)=zthmin1
      zBa(ISRCH)=zBmin1
      call rfgbmmx1(ztha,zBa,ISRCH,1,zthBmin,zBmin)

      ztha(1)=zthmax0
      zBa(1)=zBmax0
      ztha(ISRCH)=zthmax1
      zBa(ISRCH)=zBmax1
      call rfgbmmx1(ztha,zBa,ISRCH,2,zthBmax,zBmax)

      return
      end

c-----------------------------------------------

      subroutine rfgbmmx0(ztha,zBa,isiz,imin,imax,immx)
c
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER immx,isiz,i
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zdth,zbmin,zbmax,zth,rfgbinv
C============
      REAL*8 ztha(isiz)    ! theta values
      REAL*8 zBa(isiz)     ! B values
c
      integer imin,imax  ! min and max indices (output)
c
c  the extreme values of ztha and zBa are supplied; fill in the interior
c  points and find the min and max values.
c
c  the min and max values are expected to be in the interior region, or
c  else an error is reported and the run halted.
c
c-----------------
c
      zdth=ztha(isiz)-ztha(1)
c
      if(zBa(1).lt.zBa(isiz)) then
         imin=1
         zbmin=zba(1)
         imax=isiz
         zbmax=zba(isiz)
      else
         imin=isiz
         zbmin=zba(isiz)
         imax=1
         zbmax=zba(1)
      endif
c
      do i=2,isiz-1
         zth=ztha(1) + AREAL(i-1)*zdth/AREAL(isiz-1)
         ztha(i)=zth
         zba(i)=rfgbinv(zth)
         if(zba(i).lt.zbmin) then
            imin=i
            zbmin=zba(i)
         endif
         if(zba(i).gt.zbmax) then
            imax=i
            zbmax=zba(i)
         endif
      enddo
c
      if((imin.eq.1).or.(imin.eq.isiz)) then
         if(immx.ne.2) then
            call errmsg_exit('rfgbmmx error -- min B at boundary')
         endif
      endif
c
      if((imax.eq.1).or.(imax.eq.isiz)) then
         if(immx.ne.1) then
            call errmsg_exit('rfgbmmx error -- max B at boundary')
         endif
      endif
c
      return
      end

c-----------------------------------------------------

      subroutine rfgbmmx1(ztha,zBa,isiz,immx,zthout,zBout)

C============
C idecl:  explicitize implicit INTEGER declarations:
!
      IMPLICIT NONE
      REAL*8, PARAMETER :: ZERO = 0.0d0
      REAL*8, PARAMETER :: HALF = 0.5d0
      REAL*8, PARAMETER :: TWO = 2.0d0
!
      INTEGER immx,isiz,iters,ii,imin,imax,indx
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zthout,zbout,zth1,zth2,zb1,zb2,zb0,zth0,zdth,zdb1,zdb2,
     > zdenom,zdelth,rfgbinv
C============
      REAL*8 ztha(isiz),zBa(isiz)   ! theta and B(theta) values
c
c  extrema points ztha(1),zBa(1), ztha(isiz),zBa(isiz) ... are input
c
c  ztha and zBa arrays are used as workspace, and overwritten.
c
c  immx -- input -- =1 to look for minimum; =2 to look for maximum.
c  zthout -- output -- theta of min(max) B
c  zBout --  output -- the min(max) B itself
c
c------------------------
c
      data iters/2/
c
      zth1=ztha(1)
      zth2=ztha(isiz)
      zB1=zBa(1)
      zB2=zBa(isiz)
c
      do ii=1,iters
         ztha(1)=zth1
         zBa(1)=zB1
         ztha(isiz)=zth2
         zBa(isiz)=zB2

         call rfgbmmx0(ztha,zBa,isiz,imin,imax,immx)

         if(immx.eq.1) then
            indx=imin
         else
            indx=imax
         endif

         zth1=ztha(indx-1)
         zth2=ztha(indx+1)
         zB1=zBa(indx-1)
         zB2=zBa(indx+1)
      enddo
c
c  and evaluate final output point -- zero slope based on
c  quadratic fit to 3 pts nearest min- or maxima
c
      zB0=zBa(indx)
      zth0=ztha(indx)
      zdth=HALF*(ztha(indx+1)-ztha(indx-1))
c
      zdB1=zB1-zB0
      zdB2=zB2-zB0
      zdenom=TWO*(zdB1+zdB2)
      if(zdenom.eq.ZERO) then
         zdelth=ZERO
      else
         zdelth=zdth*(zdB1-zdB2)/zdenom
      endif
c
      zthout=zth0+zdelth
      zBout =rfgbinv(zthout)

      return
      end
