module volume_rings

      implicit none
      real, allocatable, save:: cosav(:,:),sinav(:,:)
      real, allocatable, save:: dv(:,:),vmean(:,:),span(:,:)
      real, allocatable, save:: xmid(:,:)
      real, allocatable, save:: cos_int(:),sin_int(:)

      logical, save :: printit


      end module volume_rings


!> therest.F :
!!
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> bzz(...)

!! Convert BI,BJ,BK (defined on cell-faces in grid coordinates) to 
!! BX,BY,BZ (cell-centered in carteisan coordinates)
!!
      SUBROUTINE BZZ(BX,BY,BZ,BI,BJ,BK)
! ... Local variables ..................................................
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
#ifdef SPHERE
#include "meter.inc"
#endif
C
      COMMON /SCRACH/
     $     FDUM(NI,NJ,NK),XDUM(NIP1,NJP1,NKP1),
     $     YDUM(NIP1,NJP1,NKP1),ZDUM(NIP1,NJP1,NKP1),
     $     bxsum(nk),bysum(nk),bzsum(nk),
     $     bx_out(nk),by_out(nk),bz_out(nk)
C
      DIMENSION XAVERI(NIP1,NJ,NK),XAVERJ(NI,NJP1,NK),XAVERK(NI,NJ,NKP1)
      DIMENSION YAVERI(NIP1,NJ,NK),YAVERJ(NI,NJP1,NK),YAVERK(NI,NJ,NKP1)
      DIMENSION ZAVERI(NIP1,NJ,NK),ZAVERJ(NI,NJP1,NK),ZAVERK(NI,NJ,NKP1)
      EQUIVALENCE (XAVERI,XDUM),(XAVERJ,XDUM),(XAVERK,XDUM)
      EQUIVALENCE (YAVERI,YDUM),(YAVERJ,YDUM),(YAVERK,YDUM)
      EQUIVALENCE (ZAVERI,ZDUM),(ZAVERJ,ZDUM),(ZAVERK,ZDUM)
      DIMENSION vol_out(nk), volsum(nk)
! ... Parameter variables ..............................................
      DIMENSION BX(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
      DIMENSION BY(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
      DIMENSION BZ(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
      DIMENSION BI(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
      DIMENSION BJ(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
      DIMENSION BK(ILOW:IHIGH,JLOW:JHIGH,KLOW:KHIGH)
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in therest.F::BZZ(...)"
#endif
C
      if ( .not. MHD_PE ) return
C
      DO 200 K=1,NL_K
      DO 200 J=1,NL_J
      DO 200 I=1,NL_I
 200     FDUM(I,J,K) = BI(I,J,K)*Face_I(I,J,k)
     $        +BI(I+1,J,K)*Face_I(I+1,J,K)
C
      DO 210 K=1,NL_K
      DO 210 J=1,NL_J
      DO 210 I=1,NL_I+1
      XAVERI(I,J,K) = X(I,J,K)+X(I,J+1,K)+X(I,J,K+1)+X(I,J+1,K+1)
      YAVERI(I,J,K) = Y(I,J,K)+Y(I,J+1,K)+Y(I,J,K+1)+Y(I,J+1,K+1)
  210 ZAVERI(I,J,K) = Z(I,J,K)+Z(I,J+1,K)+Z(I,J,K+1)+Z(I,J+1,K+1)
C
      DO 220 K=1,NL_K
      DO 220 J=1,NL_J
      DO 220 I=1,NL_I
      BX(I,J,K) = FDUM(I,J,K)*(XAVERI(I+1,J,K)-XAVERI(I,J,K))
      BY(I,J,K) = FDUM(I,J,K)*(YAVERI(I+1,J,K)-YAVERI(I,J,K))
  220 BZ(I,J,K) = FDUM(I,J,K)*(ZAVERI(I+1,J,K)-ZAVERI(I,J,K))
C
C
        DO 300 K=1,NL_K
        DO 300 J=1,NL_J
        DO 300 I=1,NL_I
  300   FDUM(I,J,K)=BJ(I,J,K)*Face_J(I,J,K)+BJ(I,J+1,K)*Face_J(I,J+1,K)
C
        DO 310 K=1,NL_K
        DO 310 J=1,NL_J+1
        DO 310 I=1,NL_I
        XAVERJ(I,J,K) = X(I,J,K)+X(I+1,J,K)+X(I,J,K+1)+X(I+1,J,K+1)
        YAVERJ(I,J,K) = Y(I,J,K)+Y(I+1,J,K)+Y(I,J,K+1)+Y(I+1,J,K+1)
  310   ZAVERJ(I,J,K) = Z(I,J,K)+Z(I+1,J,K)+Z(I,J,K+1)+Z(I+1,J,K+1)
C
        DO 320 K=1,NL_K
        DO 320 J=1,NL_J
        DO 320 I=1,NL_I
        BX(I,J,K)= BX(I,J,K)+FDUM(I,J,K)*(XAVERJ(I,J+1,K)-XAVERJ(I,J,K))
        BY(I,J,K)= BY(I,J,K)+FDUM(I,J,K)*(YAVERJ(I,J+1,K)-YAVERJ(I,J,K))
  320   BZ(I,J,K)= BZ(I,J,K)+FDUM(I,J,K)*(ZAVERJ(I,J+1,K)-ZAVERJ(I,J,K))
C
C
        DO 400 K=1,NL_K
        DO 400 J=1,NL_J
        DO 400 I=1,NL_I
  400   FDUM(I,J,K)=BK(I,J,K)*Face_K(I,J,K)+BK(I,J,K+1)*Face_K(I,J,K+1)
C
        DO 410 K=1,NL_K+1
        DO 410 J=1,NL_J
        DO 410 I=1,NL_I
        XAVERK(I,J,K) = X(I,J,K)+X(I+1,J,K)+X(I,J+1,K)+X(I+1,J+1,K)
        YAVERK(I,J,K) = Y(I,J,K)+Y(I+1,J,K)+Y(I,J+1,K)+Y(I+1,J+1,K)
  410   ZAVERK(I,J,K) = Z(I,J,K)+Z(I+1,J,K)+Z(I,J+1,K)+Z(I+1,J+1,K)
C
        DO 420 K=1,NL_K
        DO 420 J=1,NL_J
        DO 420 I=1,NL_I
        BX(I,J,K)= BX(I,J,K)+FDUM(I,J,K)*(XAVERK(I,J,K+1)-XAVERK(I,J,K))
        BY(I,J,K)= BY(I,J,K)+FDUM(I,J,K)*(YAVERK(I,J,K+1)-YAVERK(I,J,K))
  420   BZ(I,J,K)= BZ(I,J,K)+FDUM(I,J,K)*(ZAVERK(I,J,K+1)-ZAVERK(I,J,K))
C
C
        DO 500 K=1,NL_K
        DO 500 J=1,NL_J
        DO 500 I=1,NL_I
        FDUM(I,J,K) = 0.125/VOLUME(I,J,K)
        BX(I,J,K) = BX(I,J,K)*FDUM(I,J,K)
        BY(I,J,K) = BY(I,J,K)*FDUM(I,J,K)
  500   BZ(I,J,K) = BZ(I,J,K)*FDUM(I,J,K)
C
#ifdef SPHERE
#ifdef DO_RINGAV
C
C    AVERAGE ALONG THE X-AXIS
C
c   need to do ringav ?
      if ( j_global_off .gt. nrings+1 .and.
     $       j_global_off + nl_j .lt. nj_global - nrings-1 ) return


c
c   loop over all j's, first for cell centered quantities
c


      do j=1,nl_j
c
c   pick front side or back side
c
         jglobal = j + j_global_off
c
#ifdef RINGAV_TEST
      write (iunit,*) 'looping over j = ',j,jglobal
#endif
c
c  front side
         if  ( jglobal .le. nrings ) then
              nring_now = jglobal
c back side
         elseif ( jglobal .ge. nj_global - nrings+1) then
            nring_now = nj_global - jglobal + 1
         else
              nring_now = 0
         endif
#ifdef RINGAV_TEST
         write (iunit, *) 'ring now ', nring_now
#endif


      if ( nring_now .gt. 0 ) then
          nav = naver(nring_now)

                
C
          do i=1,nl_i
             do k=1,nl_k
                bxsum(k) = bx(i,j,k)
                bysum(k) = by(i,j,k)
                bzsum(k) = bz(i,j,k)
             enddo
           call phi_av(bxsum,bx_out,nl_k,nring_now,nav,.true.,.false.)
           call phi_av(bysum,by_out,nl_k,nring_now,nav,.true.,.false.)
           call phi_av(bzsum,bz_out,nl_k,nring_now,nav,.true.,.false.)
             do k=1,nl_k
                bx(i,j,k) = bx_out(k)
                by(i,j,k) = by_out(k)
                bz(i,j,k) = bz_out(k)
             enddo
          enddo
       endif
       enddo
c
#endif
#endif
C

C
! ... End SUBROUTINE BZZ(BX,BY,BZ,BI,BJ,BK) .............................
        RETURN
        END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> cross(x1,x2)
!! Take Cross product of two dimensional vectors.
!! @return cross = (x1) x (x2)
!!
        function cross(x1,x2)
! ... Parameter variables ..............................................        
        dimension x1(2)
!> x1 - first 2d vector to take cross product with.       
        dimension x2(2)
!> x2 - second 2d vector to take cross product with.       
! ... Begin ............................................................
        cross = x1(1)*x2(2) - x2(1)*x1(2)
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#ifdef SPHERE
!> ringav
!!
!!
        subroutine ringav
! ... Local variables ..................................................
        use volume_rings
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
#include "run-time.inc"
#include "help.inc"
#include "meter.inc"
#include "b_dbl.inc"
C
        common /scrach/ DUMMY(NIP1,NK,12)
c
      DIMENSION dsum(nk),rhosum(nk),vxsum(nk),vysum(nk),
     $          vzsum(nk),csum(nk),bisum(nk),
     $          bjsum(nk),volsum(nk)
      DIMENSION rho_out(nk),vx_out(nk),vy_out(nk),vz_out(nk),c_out(nk),
     $     bi_out(nk),bj_out(nk),vol_out(nk),
     $     ei_star(ni,nkp1),ej_star(nip1,nkp1)
c
#ifdef RINGAV_TEST
      dimension voltest(nk)
#endif
c
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in therest.F::ringav(...)"
#endif
        nl_Ip1 = nl_i + 1
        nl_jp1 = nl_j + 1
        nl_kp1 = nl_k + 1
#ifdef RINGAV_TEST
        iunit = 20+mype
        write (iunit,*) 'start of ringav on processor  ', mype
        write (iunit,*) 'global dims',ni_global,nj_global,nk_global
        write (iunit,1887) i_global_off,nl_i,ilow,ihigh,iactive
        write (iunit,1887) j_global_off,nl_j,jlow,jhigh,jactive
        write (iunit,1887) k_global_off,nl_k,klow,khigh,kactive
        flush(iunit)
 1887   format(10i7)
#endif
C     
c
C    AVERAGE ALONG THE X-AXIS
C

c   need to do ringav ?
      if ( j_global_off .gt. nrings+1 .and.
     $       j_global_off + nl_j .lt. nj_global - nrings-1 ) return


c
c   loop over all j's, first for cell centered quantities
c


      do j=1,nl_j
c
c   pick front side or back side
c
         jglobal = j + j_global_off
c
#ifdef RINGAV_TEST
      write (iunit,*) 'looping over j = ',j,jglobal
#endif
c
c  front side
         if  ( jglobal .le. nrings ) then
              nring_now = jglobal
c back side
         elseif ( jglobal .ge. nj_global - nrings+1) then
            nring_now = nj_global - jglobal + 1
         else
              nring_now = 0
         endif
#ifdef RINGAV_TEST
         write (iunit, *) 'nring now ', nring_now
#endif


      if ( nring_now .gt. 0 ) then
          nav = naver(nring_now)


          do i=1,nl_i
             do k=1,nl_k
                rhosum(k) = rho_g(i,j,k)
                vxsum(k) = vx_g(i,j,k)
                vysum(k) = vy_g(i,j,k)
                vzsum(k) = vz_g(i,j,k)
                csum(k) = rhosum(k)*c_g(i,j,k)**2
             enddo
           call phi_av(rhosum,rho_out,nl_k,nring_now,nav,
     $            .true.,.false.)
           call phi_av(vxsum,vx_out,nl_k,nring_now,nav,.true.,.false.)
           call phi_av(vysum,vy_out,nl_k,nring_now,nav,.true.,.false.)
           call phi_av(vzsum,vz_out,nl_k,nring_now,nav,.true.,.false.)
           call phi_av(csum,c_out,nl_k,nring_now,nav,.true.,.false.)
#ifdef RINGAV_TEST 
             if (i .eq. 8 ) then
                do kz = 1,nk
                   voltest(kz) = volume(i,j,kz)
                enddo
             iunit = 20 +mype
             write(iunit,*)'new ring ', mype,
     $            j,j+jglobal_off,nring_now,nav
             call print_avs('rho',rhosum,rho_out,voltest,nl_k,nav,iunit)
             call print_avs('vx',vxsum,vx_out,voltest,nl_k,nav,iunit)
             call print_avs('vy',vysum,vy_out,voltest,nl_k,nav,iunit)
             call print_avs('vz',vzsum,vz_out,voltest,nl_k,nav,iunit)
             call print_avs('c',csum,c_out,voltest,nl_k,nav,iunit)
             endif
             flush(iunit)
#endif
             do k=1,nl_k
                rho_g(i,j,k) = rho_out(k)
                vx_g(i,j,k) = vx_out(k)
                vy_g(i,j,k) = vy_out(k)
                vz_g(i,j,k) = vz_out(k)
                c_g(i,j,k) = sqrt(max(1.,c_out(k)/rho_out(k)))
             enddo
          enddo
c
c

          do i=1,nl_ip1
c
c             
             do k=1,nl_kp1
                ej_star(i,k) = 0.
             enddo
             do k=1,nl_k 
                bisum(k) = bi_g(i,j,k)
             enddo
             nav = naver(nring_now)
           call phi_av(bisum,bi_out,nl_k,nring_now,nav,.false.,.true.)
#ifdef RINGAV_TEST
             if ( i .eq. 8 ) then
             call print_avs('bi',bisum,bi_out,voltest,nl_k,nav,iunit)
             endif
#endif   
             do k=1,nl_k,nav
                do kl = 1, nav-1
                  ej_star(i,k+kl)=(-bi_out(kl+k-1) + bi_g(i,j,kl+k-1))
     $                  *face_i(i,j,kl+k-1)
     $                +  ej_star(i,k+kl-1)
                enddo
             enddo
c
c
c  construction allows ringav in standalone mode
c           dt > 0 => normal loop
c           dt = 0 => standalone
c           flux_fix must be run afterward in both cases to get mag fields right
c
           if ( dt .gt. 0.0) then
              do k=1,nl_kp1
                 ej_g(i,j,k) = ej_g(i,j,k)  +
     $                ej_star(i,k)/dt
              enddo
           else
              do k=1,nl_kp1
                 ej_g(i,j,k) = ej_star(i,k)
              enddo
           endif

           enddo
c
        endif  !! if check for nring_now in bounds
c     
        enddo  !! closing the original j loop
c
c


      do j=1,nl_jp1
c
c   pick front side or back side
c
         jglobal = j + j_global_off
c
c  front side
         if  ( jglobal .le. nrings+1 .and. jglobal .gt. 1) then
              nring_now = jglobal-1
c back side
         elseif ( jglobal .ge. nj_global - nrings+1 .and.
     $        jglobal .le. nj_global) then
            nring_now = nj_global - jglobal+1
         else
              nring_now = 0
         endif
         jb = j
#ifdef RINGAV_TEST
         write (iunit, *) 'nring now ', nring_now,jb
#endif
c
c
      if ( nring_now .gt. 0 ) then
           nav = naver(nring_now)
           do i=1,nl_i
c
c
              do k=1,nl_k
                 bjsum(k) = bj_g(i,jb,k)
                 ei_star(i,k) = 0.
              enddo
                 ei_star(i,nl_kp1) = 0.
              nav = naver(nring_now)
            call phi_av(bjsum,bj_out,nl_k,nring_now,nav,.false.,.true.)
#ifdef RINGAV_TEST 
             if (i .eq. 8 ) then
                do kz = 1,nk
                   voltest(kz) = face_j(i,jb,kz)
                enddo
             iunit = 20 +mype
             call print_avs('bj',bjsum,bj_out,voltest,nl_k,nav,iunit)
             endif
             flush(iunit)
#endif
              do k=1,nl_k,nav
                do kl=1,nav-1
                   ei_star(i,k+kl) = (bj_out(kl+k-1)-bj_g(i,jb,kl+k-1))*
     $                  face_j(i,jb,kl+k-1) 
     $                  + ei_star(i,k+kl-1)
                enddo
              enddo
            enddo
c
c
c
c  construction allows ringav in standalone mode
c           dt > 0 => normal loop
c           dt = 0 => standalone
c           flux_fix must be run afterward in both cases to get mag fields right
c
          if ( dt .gt. 0.0 ) then  
            do i=1,nl_i
               do k=1,nl_kp1
                  ei_g(i,jb,k) = ei_g(i,jb,k)  +
     $                 ei_star(i,k)/dt
               enddo
            enddo
          else  
            do i=1,nl_i
               do k=1,nl_kp1
                  ei_g(i,jb,k) =  ei_star(i,k)
               enddo
            enddo
          endif
      endif

      enddo

!     
! ... End subroutine ringav .............................................
!
        return
        end
#endif
!      end of "SPHERE" ifdef
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> ave2flux
!!
        subroutine ave2flux
! ... Local variables ..................................................
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in therest.F::ave2flux()"
#endif
C
      if ( .not. MHD_PE ) return
C 
        DO 340 K=1,NL_K
        DO 340 J=1,NL_J
        DO 340 I=1,NL_I+1
  340   BI_g(I,J,K) = BI_g(I,J,K)*Face_I(I,J,K)
C
        DO 350 K=1,NL_K
        DO 350 J=1,NL_J+1
        DO 350 I=1,NL_I
  350   BJ_g(I,J,K) = BJ_g(I,J,K)*Face_J(I,J,K)
C
        DO 360 K=1,NL_K+1
        DO 360 J=1,NL_J
        DO 360 I=1,NL_I
  360   BK_g(I,J,K) = BK_g(I,J,K)*Face_K(I,J,K)
!
! ... End subroutine ave2flux ...........................................
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> flux2ave
!!
!!
        subroutine flux2ave
! ... Local variables ..................................................
#include "global_dims.inc"
#include "mhd_P++_pointers.inc"
#include "fortran_P++_arrays.inc"
#include "param.inc"
! ... Begin ............................................................
#ifdef DEBUG_MODE_ON
      write(*,*) "DEBUG:  in therest.F::flux2ave()"
#endif
C
      if ( .not. MHD_PE ) return
C 
        DO 340 K=1,NL_K
        DO 340 J=1,NL_J
        DO 340 I=1,NL_I+1
  340   BI_g(I,J,K) = BI_g(I,J,K)/Face_I(I,J,K)
C
        DO 350 K=1,NL_K
        DO 350 J=1,NL_J+1
        DO 350 I=1,NL_I
  350   BJ_g(I,J,K) = BJ_g(I,J,K)/Face_J(I,J,K)
C
        DO 360 K=1,NL_K+1
        DO 360 J=1,NL_J
        DO 360 I=1,NL_I
  360   BK_g(I,J,K) = BK_g(I,J,K)/Face_K(I,J,K)
!
! ... End subroutine flux2ave ...........................................
        return
        end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c
      subroutine phi_av(FIN, FOUT, N, NR,NAV,use_slope,fourier_reduce)
c
c   averages over phi direction using MUSCL decomposition
c   fin is cell integral \int f dV input
c   fout is the averaged value
c   n dimension in the phi direction
c   nav number of cells to average over
c
c
c
      use volume_rings

      dimension fin(n), fout(n)
      logical fourier_reduce, use_slope
      real, allocatable:: fstar(:), fav(:)
      real, allocatable:: fout_av(:), fs_out(:), fs_diff(:)
      real, allocatable:: slope(:),slope2(:),sl_sign(:)
      real, allocatable:: fourier(:)
c

      twopi = 8*atan(1.0)
      allocate(fourier(n))
      allocate(fstar(n))
      allocate(fs_out(n))
#ifdef PHIAV_TEST
      allocate(FS_diff(n))
      allocate(fout_av(n))
#endif
      ncells = n/nav
c ncells = number of averaged regions
      allocate(fav(ncells))
      allocate(slope(ncells))
      allocate(slope2(0:ncells))
      allocate(sl_sign(0:ncells))
c
c total mass with in the ring 
      fsum = 0.
      do i=1,n
         fsum = fsum + fin(i)
      enddo
c
      piinv = 1./(4.*atan(1.0))
c
      if ( fourier_reduce) then
         cos_coef = 0.
         sin_coef =0.
         do i=1,n
            cos_coef = cos_coef + fin(i)*cos_int(I)
            sin_coef = sin_coef + fin(i)*sin_int(i)
          enddo
          cos_coef = piinv*cos_coef
          sin_coef = piinv*sin_coef
c
          do i = 1,n
             fourier(i) =
     $            (cos_coef*cos_int(I)+sin_coef*sin_int(i))/dv(i,nr)
             fin(i) = fin(i) - fourier(i)
          enddo
       endif



      do i=1,ncells
         fav(i) = 0.
         do j=1,nav
            fav(i) = fav(i) + fin((i-1)*nav+j)*dv((i-1)*nav+j,nr)
         enddo
         fav(i) = fav(i)/vmean(i,nr)
      enddo
      do i=1,ncells
         do j=1,nav
            fout((i-1)*nav+j) = fav(i)
         enddo
      enddo
c
      if (use_slope) then
      slope2(0) = fav(1)-fav(ncells)
      slope2(ncells) = fav(1)-fav(ncells)
      do i=1,ncells-1
         slope2(i) = fav(i+1)-fav(i)
      enddo
      do i = 0,ncells
         sl_sign(i) = sign(1.0,slope2(i))
         if (slope2(i) .eq. 0.0 ) sl_sign(i) = 0.0 
      enddo

      slope(1) = 0.5*(fav(2)-fav(ncells))
      slope(ncells) = 0.5*(fav(1)-fav(ncells-1))
      do i=2,ncells-1
         slope(i) = 0.5*(fav(i+1) - fav(i-1))
      enddo
c
      do i=1,ncells
         if ( sl_sign(i-1)*sl_sign(i) .le. 0.0 ) then
            slope(i) = 0.
         else
            slope(i) = sign(min(abs(slope(i)),
     $           abs(2.*slope2(i-1)),
     $           abs(2.*slope2(i))),
     $           slope(i))
         endif
         slope(i) = slope(i)/float(nav)
      enddo
      do i=1,ncells
         do j=1,nav
            fout((i-1)*nav+j) = fout((i-1)*nav+j)
     $        +   slope(i)*(j-xmid(i,nr))
         enddo
      enddo
      endif
c
c
      if (fourier_reduce) then
         do i =1,n
            fout(i) = fout(i) + fourier(i)
            fin(i) = fin(i) + fourier(i)
         enddo
      endif
c
c
#ifdef PHIAV_TEST
c  testing to see whether decomposition is working
c  print input and output values
c
      do i=1,n
         fs_diff(i) = fstar(i) - fs_out(i)
      enddo

c
      do i=1,n
      write (6,901) i, fin(i),fout(i),fstar(i),fs_out(i),fs_diff(I)
 901  format (i5,6(1pe12.4))
      enddo
c
c  compute the averaged region vlues for input and output and
c  print the comparison
c
      do i=1,ncells
         fav(i) = 0.
         fout_av(i) = 0.
         do j=1,nav
            fout_av(i) = fout_av(i) + fout((i-1)*nav+j)
            fav(i) = fav(i) + fin((i-1)*nav+j)
         enddo
      enddo

      do j=1,ncells
      write (6,902) j, fav(j),fout_av(j),fav(j)-fout_av(j)
 902  format (i4, 3(1pe12.4))
      enddo
      deallocate(fs_diff)
      deallocate(fout_av)
#endif
      deallocate(fourier)
      deallocate(sl_sign)
      deallocate(slope2)
      deallocate(slope)
      deallocate(fstar)
      deallocate(fs_out)
      deallocate(fav)

c
      return
      end
c
c
#ifdef RINGAV_TEST
      subroutine print_avs(title,fin,fout,vol,nk,nav,io)
c
      character*(*) title
      real fin(nk),fout(nk),vol(nk)
      real, allocatable:: fav(:),fout_av(:),fin_tot(:),fout_tot(:)
c
      ncells = nk/nav
c
      allocate(fav(ncells))
      allocate(fout_av(ncells))
      allocate(fin_tot(nk))
      allocate(fout_tot(nk))
c  compute the averaged region vlues for input and output and
c  print the comparison
c
      do i=1,nk
         fin_tot(i) = fin(i)*vol(i)
         fout_tot(i) = fout(i)*vol(i)
      enddo
      do i=1,ncells
         fav(i) = 0.
         fout_av(i) = 0.
         do j=1,nav
            fout_av(i) = fout_av(i) + fout_tot((i-1)*nav+j)
            fav(i) = fav(i) + fin_tot((i-1)*nav+j)
         enddo
      enddo
      write (io,901) title
 901  format(A)
      do j=1,ncells
      write (io,902) j, fav(j),fout_av(j),fav(j)-fout_av(j)
 902  format (i4, 3(1pe12.4))
      enddo    
      deallocate(fav)
      deallocate(fout_av)
      deallocate(fin_tot)
      deallocate(fout_tot)
      flush(io)
c
      return
      end
#endif
c

      subroutine ring_setup(fin,n,nrings,naver)
      use volume_rings
#include "global_dims.inc"

      real fin(n)
      integer naver(n)
c
      if (nrings .eq. 0 ) return
      navmin = naver(nrings)
      ncell_max = n/navmin
c
         allocate(cos_int(n))
         allocate(sin_int(n))
         allocate(cosav(n,nrings))
         allocate(sinav(n,nrings))
         allocate(dv(n,nrings))
         allocate(vmean(ncell_max,nrings))
         allocate(span(ncell_max,nrings))
         allocate(xmid(ncell_max,nrings))
c
c
         twopi = 8.0*atan(1.0)
         voltot = 0.
         do i=1,n
            voltot = voltot + fin(i)
         enddo
c
         philow = 0.
         do i=1,n
            phihigh = philow + twopi*fin(i)/voltot
            cos_int(i) = sin(phihigh)-sin(philow)
            sin_int(i) = cos(philow)-cos(phihigh)
            philow = phihigh
         enddo
c
c

c
c
         do nr = 1, nrings
c
         nav = naver(nr)
         ncells = n/nav
c
         do i=1,n
            dv(i,nr) = fin(i)*twopi/voltot
         enddo
         do j=1,ncells
            vmean(j,nr) = 0.
            visum = 0.
            do i=1,nav
               vmean(j,nr) = vmean(j,nr) + dv((j-1)*nav+i,nr)
               visum = visum + dv((j-1)*nav+i,nr)*float(i)
            enddo
               xmid(j,nr) = visum/vmean(j,nr)
               sum_test = 0.
               do i=1,nav
                  sum_test = sum_test + dv((j-1)*nav+i,nr)*
     $                 (float(i)-xmid(j,nr))
               enddo
         enddo
c
c
c
c  do bj stuff
c
         do j = 1,ncells
            span(j,nr) = 1./(2.*sin(vmean(j,nr)*0.5))
            phi0 = vmean(j,nr)*0.5
               cosav((j-1)*nav+1,nr) = (sin(-phi0+dv((j-1)*nav+1,nr))-
     $           sin(-phi0))
               sinav((j-1)*nav+1,nr) = -cos(-phi0+dv((j-1)*nav+1,nr))+
     $           cos(-phi0)
               philow = -phi0 + dv((j-1)*nav+1,nr)
            do i=2,nav
               cosav((j-1)*nav+i,nr) = (sin(philow+dv((j-1)*nav+i,nr))-
     $           sin(philow))
               sinav((j-1)*nav+i,nr) = -cos(philow+dv((j-1)*nav+i,nr))+
     $           cos(philow)
               philow = philow+dv((j-1)*nav+i,nr)
            enddo
         enddo
c       
c         
         enddo
         return
         end