#include "fpreproc/f77_dcomplx.h"
C
c  dmc 24 Oct 1997 -- completing multi-species generalization.
c
c  ** rfxqlo_psical is called outside of the species loop:  therefore,
c     dependence of quantities on rZbeam or Ab are not allowed, so:
c
c  ** taking rZbeam/Ab factor out of definition of omppofr(...) array,
c     and factoring this in where omppofr is used (thupdown subroutine
c     under rfoper).
c
c  ** taking rZbeam factor out of PSI
c     and putting it back where these factors are used
c     (lossco,rfoper)
c
c  dmc 21 Jan 1997 at TEXTOR --
c   starting geometry generalization
c   still using B=k/R; updown symmetric plasma can now have indentation
c   (as in PBXM) but updown *asymmetric* plasma should have only one
c   magnetic well -- i.e. only one Bmin and one Bmax --
c   on all surfaces -- else theta(B) root finder (r8_zeroin) will fail
c   because of singularities dB/dtheta = 0.
c   ...more work will be needed to support more general field topologies
c
c   use rfgbmmx subroutine to find theta(Bmin),Bmin, theta(Bmax),Bmax
c
c---------------------
c  dmc June 2003 -- further adjustments:
c    (a) remove remnants of circular geometry approximation; use
c        x=sqrt(phi/philim) instead of r in calculations.
c    (b) assume existance of unique Bmin(Rmax) and Bmax(Rmin);
c        B(R) representation will be "simplified" to conform to
c        this as needed.
c---------------------
c
C	THIS PROGRAM CALCULATES THE POLOIDAL FLUX PROFILE
C	FOR THE BOUNCE-AVERAGED FOKKER-PLANCK CODE
C
C	IT FILLS AN ARRAY WITH PSI(R), AND A SECOND
C	ARRAY WITH R(PSI) FOR USE IN THE FINITE BANANA
C	WIDTH INJECTION, DETECTION, AND LOSS-CONE SECTIONS
C
C	PSI(1) IS PSI @ X=0 (= 0), PSI(numpsi) IS PSI @ X=1
C
C	IT USES q(r) as specified in the array qzc.
c
c--Particle orbits in a tokamak can be calculated from
c	conservation of canonical angular momentum:
c	Pphi = m R vparallel - psiec = constant
c
c	psiec = -Aphi R e/c
c
c	d(psiec)/dr = (e/c)*R*Bp
C
C	ROB GOLDSTON OCT. 4, 1981
C
C---------------------------------------
C
      SUBROUTINE RFXQLO_PSICAL
C
C	COMMON BLOCKS
C
      use rfxqlo_private

C============
C idecl:  explicitize implicit INTEGER declarations:
!
      IMPLICIT NONE
      REAL*8, PARAMETER :: ZERO = 0.0d0
      REAL*8, PARAMETER :: HALF = 0.5d0
      REAL*8, PARAMETER :: C100TH = 0.01d0
      REAL*8, PARAMETER :: ONE = 1.0d0
      REAL*8, PARAMETER :: TWO = 2.0d0
      REAL*8, PARAMETER :: CEPS6 = 1.0d-6
      REAL*8, PARAMETER :: CBIG6 = 1.0d6
      REAL*8, PARAMETER :: CBIG8 = 1.0d8
!
      INTEGER i,iz,ir,j,isides,is,jb,imaxit
      integer ierr
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zzeps,eovc,zpi,zb,zxloc,zqloc,rfintrp,dphidx,
     > zthbmin,zbmin,zthbmax,zbmax,rfgbinv,zrmajmax,
     > zybmin,zybmax,
     > zdthstar,zdthact,zth1,zrdum,zydum,zj0,zb0,zdth,zjp,zdrp,zjm,
     > zdrm,zthstar,zth0,r8_zeroin
      REAL*8 zrmin,zrmax,zzmin,zzmax,zdum
C============
      external rfgbinv ! calling root finder with this routine will
      REAL*8 zBrequ ! find the theta giving this zBrequ
      REAL*8 zXi    ! on this flux surface (normalized rho/rho_lim)
      REAL*8 zRB    ! R*Bt for rfgbcom
      common /rfgbcom/ zxi, zBrequ, zRB
      REAL*8 zJac(2,2)	! jacobian calculated by MOMRY
c
      REAL*8 zxa(2*nx_qlo)
c
c  dmc 22 Feb 1993
c  replacing IMSL root finder ZBREN with call to R8_ZEROIN
c
      data zzeps/1.2d-7/
c
C
c--eovc = (e/c) = electron charge over the speed of light (cgs units):
      DATA EOVC/1.6021d-20/
C---------------------------------------
C  executable code
C
      zpi=abs(acos(-ONE))
      Zb=1.0		! use this instead of rZbeam dmc Oct 97
c
c [dqdrpro -- no longer needed, no longer calculated -- dmc June 2003]
c
C calculate Psi(x) from
c
c d(Psi)/dx = (ze/c 2 pi) (1/q) d(Phi)/dx
c
c   with
c
c where Phi is the toroidal flux and Psi is our definition
c of poloidal flux. (dPsi/dr = (1/q) dPhi/dr is the usual definition
c but we have an extra factor in our Psi.)
C
C note that Psi has been rescaled to be ze/c/(2 pi) times the raw
C poloidal flux.

      delxPsi=ONE/(numpsi-ONE)
      PSI(1)=ZERO
      DO 100 I=1,numpsi-1
         zxloc=(I-HALF)/(numpsi-ONE)
         zqloc=rfintrp(zxloc,qzc,x_qlo,nx_qlo)

c The following expression for d(Phi)/dr assumes that the toroidal flux
c scales as Phi(r) = (r/a)**2, which is true in cases where the radial
c grid was chosen to be uniform in square-root-of-the-toroidal-flux.
c It would be useful to upgrade this eventually to use the actual Phi(r)
c profile to allow for more general cases...
         dPhidx=CBIG8*tflux*2*zxloc

         PSI(I+1)=PSI(I)
     1     +Zb*eovc/(2*zpi)/zqloc*dPhidx*delxPsi

100   CONTINUE

c dmc -- moved the following code up out of the 400 loop, and
c  generalized for updown asymmetry.

      do i=1,nx_qlo
         do iz=0,1
            ir=2*i-1+iz
            if(iz.eq.0) then
               zxi=x_qlo(i)
               zRB=gzc(i)*CBIG6         ! to Gauss*cm
            else if(iz.eq.1) then
               zxi=xb_qlo(i)
               zRB=gzb(i)*CBIG6         ! to Gauss*cm
            endif
            zxa(ir)=zxi
            zdum=ZERO
c  need Gauss, not  T; need Th(Rmin),Th(Rmax) not Th(Bmin),Th(Bmax)...
cxx            call eq_glimb1(zxi,zdum,zBmin,zthBmin,zBmax,zthBmax,ierr)
cxx            if(ierr.ne.0)
cxx     >         call errmsg_exit('??fppmod(psical)->eq_glimb1 ierr set.')
cxx            call eq_glimrz(zxi,zRmin,zRmax,zZmin,zZmax,ierr)
cxx            if(ierr.ne.0)
cxx     >         call errmsg_exit('??fppmod(psical)->eq_glimrz ierr set.')
            zxi=min(0.9999999d0,zxi)
c  use rfgbmmx to find the min/max B and corresponding theta values
            call rfgbmmx(zthBmin,zBmin,zthBmax,zBmax)
            if(rfx_grid%symflag) then
c  set theta bounds to 0 and pi (the normal case for ud symmetry) if
c  we are close...
               if(abs(zthbmin).le.C100TH) then
                  zthbmin=ZERO
                  zbmin=rfgbinv(zthbmin)
               endif
               if((abs(zthbmax)-zpi).le.C100TH) then
                  zthbmax=zpi
                  zbmax=rfgbinv(zthbmax)
               endif
            endif
            Bminofr(ir)=zBmin
            ThBmin(ir)=zthBmin
            call rfmomry1(zxi,zthBmin,0,RatBmin(ir),YatBmin(ir),zjac)
            if(rfx_grid%symflag) then
               if(zthbmin.eq.ZERO) yatBmin(ir)=ZERO
               Bmaxofr(ir)=zBmax
               ThBmaxP(ir)=abs(zthBmax)
               ThBmaxM(ir)=-ThBmaxP(ir)
               call rfmomry1(zxi,abs(zthBmax),0,
     1            RatBmaxP(ir),YatBmaxP(ir),zjac)
               if(zthBmax.eq.zpi) YatBmaxP(ir)=ZERO
               RatBmaxM(ir)=RatBmaxP(ir)
               YatBmaxM(ir)=-YatBmaxP(ir) ! allow PBX geometry
            else
               Bmaxofr(ir)=zBmax
               if(zthBmax.gt.zthBmin) then
                  ThBmaxP(ir)=zthBmax
                  ThBmaxM(ir)=zthBmax-TWO*zpi
               else
                  ThBmaxM(ir)=zthBmax
                  ThBmaxP(ir)=zthBmax+TWO*zpi
               endif
               call rfmomry1(zxi,thBmaxP(ir),0,
     1            RatBmaxP(ir),YatBmaxP(ir),zjac)
               RatBmaxM(ir)=RatBmaxP(ir)
               YatBmaxM(ir)=YatBmaxP(ir) ! UDA:  assume reach same max
            endif
            Bavgofr(ir)=(Bmaxofr(ir)+Bminofr(ir))/TWO
            Bdelofr(ir)=(Bmaxofr(ir)-Bminofr(ir))/TWO
         enddo
      enddo
c
c calculate arrays to do a quick inverse map from (B,r) to theta
c (poloidal angle), or from (B,r) to Y (vertical position).
c This is used if you want to know where a particle's banana tip
c is, or where its resonance is.
c
c mod dmc:  support updown asymmetry.  The mappings can be different
c  btw the lower and upper halfplanes
c
c ir index runs over zone centers and zone boundaries:
c dmc:  I evaluate both up & down, even in the symmetric case;
c   this is useful for testing and makes coding easier downstream.
c   cost is small.

      isides=2

      do is=1,isides
         do i=1,nx_qlo
            do iz=0,1
               ir=2*i-1+iz
               if(iz .eq. 0) then
                  zxi=x_qlo(i)
                  zrb=gzc(i)*CBIG6         ! to Gauss*cm
               else if(iz .eq. 1) then
                  zxi=xb_qlo(i)
                  zrb=gzb(i)*CBIG6         ! to Gauss*cm
               endif
c ensure MOMRY's jacobian is evaluated just inside the boundary:
               if(zxi .ge. ONE) zxi=0.999999d0
               zthbmin=thbmin(ir)
               zybmin=yatbmin(ir)
               if(is.eq.1) then
                  zthbmax=thbmaxp(ir)
                  zybmax=yatbmaxp(ir)
               else
                  zthbmax=thbmaxm(ir)
                  zybmax=yatbmaxm(ir)
               endif

               if(is.eq.1) then
                  numthofB(ir)=maxthofb*sqrt(AREAL(i)/AREAL(nx_qlo))
                  if(numthofB(ir) .lt. 3) numthofB(ir)=3
c  dtheta(star):
                  zdthstar=zpi/(numthofB(ir)-1)
                  delthofB(ir)=zdthstar
               else
                  zdthstar=delthofB(ir)
               endif

c  moments coords dth:  caution:  this is negative for lower halfplane:
               zdthact=(zthbmax-zthbmin)/(numthofB(ir)-1)

               ThofBr(1,ir,is)=zthbmin
               ThofBr(numthofB(ir),ir,is)=zthbmax

               YofBr(1,ir,is)=zybmin
               YofBr(numthofB(ir),ir,is)=zybmax

c calculate omppofr = d**2 Omega / d l**2 at theta=0 and pi:
c   (mod dmc:  actually at theta(Bmin) and theta(Bmax) now):

               zqloc=rfintrp(zxi,qzc,x_qlo,nx_qlo)
               jb=1
               zth1=zthbmin

               call rfmomry1(zxi,zth1,1,zRdum,zYdum,zJac)
               zj0=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zB0=Zrb/zRdum
               zdth=zdthact
               call rfmomry1(zxi,zth1+zdth/TWO,1,zRdum,zYdum,zJac)
               zjp=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zdRp=zjac(1,2)/zRdum
               call rfmomry1(zxi,zth1-zdth/TWO,1,zRdum,zYdum,zJac)
               zjm=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zdRm=zjac(1,2)/zRdum
c  Z/A factor removed dmc Oct 1997:
               omppofr(jb,ir,is)=-9.579d3/zB0/zj0*
     1            (zdRp/zjp-zdRm/zjm)/zdth

               jb=numthofB(ir)
               zth1=zthbmax

               call rfmomry1(zxi,zth1,1,zRdum,zYdum,zJac)
               zj0=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zB0=Zrb/zRdum
               zdth=zdthact
               call rfmomry1(zxi,zth1+zdth/TWO,1,zRdum,zYdum,zJac)
               zjp=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zdRp=zjac(1,2)/zRdum
               call rfmomry1(zxi,zth1-zdth/TWO,1,zRdum,zYdum,zJac)
               zjm=zpi*zqloc/zxi/CBIG8/tflux*zRdum*(zjac(1,1)*zjac(2,2)-
     1            zjac(2,1)*zjac(1,2))
               zdRm=zjac(1,2)/zRdum
c  Z/A factor removed dmc Oct 1997:
               omppofr(jb,ir,is)=-9.579d3/zB0/zj0*
     1            (zdRp/zjp-zdRm/zjm)/zdth

c loop around from theta=0 to pi to find the various inverse maps:

               do 390 jb=2,numthofB(ir)-1
                  zthstar=(jb-1)*zdthstar
                  zBrequ=Bavgofr(ir)-Bdelofr(ir)*cos(zthstar)
                  zth0=min(zthbmin,zthbmax)
                  zth1=max(zthbmin,zthbmax)
                  imaxit=100
c finds theta to an accuracy of +/- 1.e-5*1.e-1:
                  zth1=r8_zeroin(zth0,zth1,rfgbinv,CEPS6,zzeps)
CIMSL	call zbren(rfgbinv,0.0,1.e-5,zth0,zth1,imaxit)
                  ThofBr(jb,ir,is)=zth1
                  call rfmomry1(zxi,zth1,1,zRdum, YofBr(jb,ir,is),zJac)

c calculate omppofr = d**2 Omega / d l**2:
                  zj0=zpi*zqloc/zxi/CBIG8/tflux*zRdum*
     1               (zjac(1,1)*zjac(2,2)-zjac(2,1)*zjac(1,2))
                  zB0=Zrb/zRdum
                  zdth=zdthact
                  call rfmomry1(zxi,zth1+zdth/TWO,1,zRdum,zYdum,zJac)
                  zjp=zpi*zqloc/zxi/CBIG8/tflux*zRdum*
     1               (zjac(1,1)*zjac(2,2)-zjac(2,1)*zjac(1,2))
                  zdRp=zjac(1,2)/zRdum
                  call rfmomry1(zxi,zth1-zdth/TWO,1,zRdum,zYdum,zJac)
                  zjm=zpi*zqloc/zxi/CBIG8/tflux*zRdum*
     1               (zjac(1,1)*zjac(2,2)-zjac(2,1)*zjac(1,2))
                  zdRm=zjac(1,2)/zRdum
c  Z/A factor removed dmc Oct 1997:
                  omppofr(jb,ir,is)=-9.579d3/zB0/zj0*
     1               (zdRp/zjp-zdRm/zjm)/zdth

 390           continue
            enddo
         enddo
      enddo

      return
      end
