      subroutine rfx_thupdown(ishe,
     >   zra,zpi,zBrequ,zthdsra,zthdsrb,zomppa,zomppb)
c
      use rfxqlo_private
c
c  dmc 24 Oct 1997:  removed rZbeam/Ab factor from omppofr array dimension
c         and returned it here (as passed argument zra)
c
c  dmc -- return the geometrical theta of the upper and lower resonances
c         and the d**2 Omega / dl**2 at the upper and lower resonances
c
c  input:  zpi, zBrequ=B of resonance; ishe is the current shell.
c
c  output:  theta angles zthdsra,zthdsrb
c           derivatives  zomppa, zomppb
c
c  we are interpolating on results previously tabulated vs. theta(star)
c  (B-related theta) previously tabulated in subroutine psical.
c
c  ** internal use only ** -- cgs units
c----------------------------------
C============
C idecl:  explicitize implicit INTEGER declarations:
!
      IMPLICIT NONE
      REAL*8, PARAMETER :: ZERO = 0.0d0
      REAL*8, PARAMETER :: ONE = 1.0d0
      REAL*8, PARAMETER :: TWO = 2.0d0
!
      INTEGER ir,isides,is,ishe
C============
C idecl:  explicitize implicit REAL declarations:
      REAL*8 zpi,zbrequ,zthdsra,zthdsrb,zomppa,zomppb,zra,zthstar,
     > zarg,rfinteq1d,z2pi
C============
      REAL*8 zth(2),zompp(2)
c
      ir=ishe*2-1
      if(rfx_grid%symflag) then
         isides=1
      else
         isides=2
      endif

      do is=1,isides
         if(zBrequ .le. Bminofr(ir)) then
            zth(is)=thBmin(ir)
            zthstar=ZERO
         else if(zBrequ .ge. Bmaxofr(ir))then
            if(is.eq.1) then
               zth(is)=thBmaxP(ir)
            else
               zth(is)=thBmaxM(ir)
            endif
            zthstar=zpi
         else
c test:
            zarg=(Bavgofr(ir)-zBrequ)/Bdelofr(ir)
            if(zarg .gt. 1.01d0.or. zarg .lt. -1.01d0)
     $         call errmsg_exit('zarg much too large')

            zthstar=acos(min(ONE,max(-ONE,zarg)))
            zth(is)=rfinteq1d(zthstar,
     $         ThofBr(1,ir,is),ZERO,delthofB(ir),numthofB(ir))
         endif
c  zra = rZbeam/Ab -- dmc Oct 1997
         zompp(is)=zra*rfinteq1d(zthstar,
     $      omppofr(1,ir,is),ZERO,delthofB(ir),numthofB(ir))
      enddo
c
      if(isides.eq.1) then
         zomppa=zompp(1)
         zomppb=zomppa
         zthdsra=zth(1)
         zthdsrb=-zth(1)
      else
         zomppa=zompp(1)
         zomppb=zompp(2)
         zthdsra=zth(1)
         zthdsrb=zth(2)
      endif
c
c  want thetas in range (0,2pi)
c
      z2pi=TWO*zpi
      if(zthdsra.lt.ZERO) then
         zthdsra=zthdsra+z2pi
      else if(zthdsra.ge.z2pi) then
         zthdsra=zthdsra-z2pi
      endif
c
      z2pi=TWO*zpi
      if(zthdsrb.lt.ZERO) then
         zthdsrb=zthdsrb+z2pi
      else if(zthdsrb.ge.z2pi) then
         zthdsrb=zthdsrb-z2pi
      endif
c
      return
      end
