!+
!NAME:
! div_3d
!PURPOSE:
! Does the divergence for a 3d vector
!CALLING SEQUENCE:
! call div_3d(ax, ay, az, x, y, z, nx, ny, nz, use_saved, $
!     &        qsphere, div_a, cxarr, cyarr, czarr)
!INPUT:
! ax, ay, az = are 3d arrays defined on (x,y,z)
! x, y, z, = the coordinates on which the arrays is defined.
! nx, ny, nz are the numbers of grid points in x, y, z
! use_saved = if set to 1, use the saved values in the deriv routines
! qsphere = if set to 1, spherical coordinates
! cxarr, cyarr, czarr are save coefficients for the derivatives
!OUTPUT:
! div_a = the divergence of a
!HISTORY:
! 22-may-2006, jmm, jimm@ssl.berkeley.edu
! 7-feb-2007, jmm, try to allocate minimal memory
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
! 22-jul-2007, jmm, Oops, bug for spherical case...
!-
      subroutine div_3d(ax, ay, az, x, y, z, nx, ny, nz, use_saved, 
     &     qsphere, div_a, cxarr, cyarr, czarr)
      implicit none
      integer :: nx, ny, nz, use_saved, qsphere
      integer :: j, k, l
      real*8 :: ax(nx, ny, nz), ay(nx, ny, nz), az(nx, ny, nz)
      real*8 :: x(nx), y(ny), z(nz), div_a(nx, ny, nz)
      real*8 :: cxarr(nx, 3), cyarr(ny, 3), czarr(nz, 3)
      real*8,allocatable :: f1(:,:,:), f0(:,:,:)
      real*8,allocatable :: sinth(:)
      allocate(f0(nx,ny,nz))
      allocate(f1(nx,ny,nz))
      If(qsphere.eq.1) then
         forall (j=1:nx,k=1:ny,l=1:nz)
            f0(j,k,l) = (x(j)**2)*ax(j,k,l)
         end forall
         call deriv_dx(x, f0, nx, ny, nz, use_saved, f1, cxarr)
         forall (j=1:nx,k=1:ny,l=1:nz)
            div_a(j,k,l) = f1(j,k,l)/(x(j)**2)
         end forall
      else
         call deriv_dx(x, ax, nx, ny, nz, use_saved, f1, cxarr)
         div_a = f1
      end if
      if(qsphere.eq.1) then
         allocate(sinth(ny))
         sinth=sin(y)
         forall (j=1:nx,k=1:ny,l=1:nz)
            f0(j,k,l) = sinth(k)*ay(j,k,l)
         end forall
         call deriv_dy(y, f0, nx, ny, nz, use_saved, f1, cyarr)
         forall (j=1:nx,k=1:ny,l=1:nz)
            div_a(j,k,l) = div_a(j,k,l)+f1(j,k,l)/(x(j)*sinth(k))
         end forall
      else
         call deriv_dy(y, ay, nx, ny, nz, use_saved, f1, cyarr)
         div_a=div_a+f1
      end if
      if(qsphere.eq.1) then
         call deriv_dz(z, az, nx, ny, nz, use_saved, f1, czarr)
         forall (j=1:nx,k=1:ny,l=1:nz)
               div_a(j,k,l) = div_a(j,k,l)+f1(j,k,l)/(x(j)*sinth(k))
         end forall
         deallocate(sinth)
      else
         call deriv_dz(z, az, nx, ny, nz, use_saved, f1, czarr)
         div_a = div_a+f1
      end if
      deallocate(f1)
      return
      end
!+
!NAME:
! curl_3d
!PURPOSE:
! Does the curl for a 3d vector
!CALLING SEQUENCE:
! call curl_3d(ax, ay, az, x, y, z, nx, ny, nz, use_ saved,
!               qsphere, curl_x, curl_y, curl_z, cxarr, cyarr, czarr)
!INPUT:
! ax, ay, az = are 3d arrays defined on (x,y,z)
! x, y, z, = the coordinates on which the arrays is defined.
! nx, ny, nz = the number of grid points in th x,y,z arrays
! use_saved = if set to 1, use the grids that have been saved
! qsphere = if set to 1, use spherical coordinates
! cxarr, cyarr, czarr are the saved grids
!OUTPUT:
! curl_x, curl_y, curl_z = the curl of a
!HISTORY:
! 18-may-2006, jmm, jimm@ssl.berkeley.edu
! 7-feb-2007, jmm, trying to save some memory
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine curl_3d(ax, ay, az, x, y, z, nx, ny, nz,
     &     use_saved, qsphere, curl_x, curl_y, curl_z, 
     &     cxarr, cyarr, czarr)
      implicit none

      integer :: nx, ny, nz, use_saved, qsphere
      integer :: j, k, l
      real*8 :: ax(nx, ny, nz), ay(nx, ny, nz), az(nx, ny, nz)
      real*8 :: x(nx), y(ny), z(nz)
      real*8 :: curl_x(nx, ny, nz), curl_y(nx, ny, nz) 
      real*8 :: curl_z(nx, ny, nz)
      real*8 :: cxarr(nx, 3), cyarr(ny, 3), czarr(nz, 3)
!all arrays that aren't passed in are allocatable
      real*8, allocatable :: f1(:,:,:)
      real*8, allocatable :: sinth(:)
!allocate memory
      allocate(f1(nx,ny,nz))
      if(qsphere.eq.1) then
         allocate(sinth(ny))
         sinth = sin(y)
!r component
         forall (j=1:nx, k=1:ny, l=1:nz)
            f1(j,k,l) = sinth(k)*az(j,k,l)
         end forall
         call deriv_dy(y, f1, nx, ny, nz, use_saved, curl_x, cyarr)
         call deriv_dz(z, ay, nx, ny, nz, use_saved, f1, czarr)
         curl_x = curl_x-f1
         forall (j=1:nx, k=1:ny, l=1:nz)
            curl_x(j,k,l) = curl_x(j,k,l)/(x(j)*sinth(k))
         end forall
!theta component
         forall (j=1:nx, k=1:ny, l=1:nz)
            f1(j,k,l) = x(j)*az(j,k,l)
         end forall
         call deriv_dx(x, f1, nx, ny, nz, use_saved, curl_y, cxarr)
         call deriv_dz(z, ax, nx, ny, nz, use_saved, f1, czarr)
         forall (j=1:nx, k=1:ny, l=1:nz)
            curl_y(j,k,l) = (f1(j,k,l)/sinth(k)-curl_y(j,k,l))/x(j)
         end forall
!phi component
         forall (j=1:nx, k=1:ny, l=1:nz)
            f1(j,k,l) = x(j)*ay(j,k,l)
         end forall
         call deriv_dx(x, f1, nx, ny, nz, use_saved, curl_z, cxarr)
         call deriv_dy(y, ax, nx, ny, nz, use_saved, f1, cyarr)
         forall (j=1:nx, k=1:ny, l=1:nz)
             curl_z(j,k,l) = (curl_z(j,k,l)-f1(j,k,l))/x(j)
         end forall
         deallocate(sinth)
      else
         call deriv_dy(y, az, nx, ny, nz, use_saved, f1, cyarr)
         curl_x = f1
         call deriv_dz(z, ay, nx, ny, nz, use_saved, f1, czarr)
         curl_x = curl_x-f1
         call deriv_dz(z, ax, nx, ny, nz, use_saved, f1, czarr)
         curl_y = f1
         call deriv_dx(x, az, nx, ny, nz, use_saved, f1, cxarr)
         curl_y = curl_y-f1
         call deriv_dx(x, ay, nx, ny, nz, use_saved, f1, cxarr)
         curl_z = f1
         call deriv_dy(y, ax, nx, ny, nz, use_saved, f1, cyarr)
         curl_z = curl_z-f1
      end if
      deallocate(f1)
      return
      end
!+
!NAME:
! grad_3d
!PURPOSE:
! Does the gradient for an array in Cartesian coordinates
!CALLING SEQUENCE:
! call grad_3d(a, x, y, z, nx, ny, nz, use_saved, qsphere, 
!    &          grad_x, grad_y, grad_z, cxarr, cyarr, czarr)
!INPUT:
! a = a 3d array
! x, y, z = the coordinates on which array is defined.
! nx, ny, nz are the numbers of grid points in x, y, z
! use_saved = if set to 1, use the saved values in the deriv routines
! qsphere = if set to 1, then spherical coordinates are used
! cxarr, cyarr, czarr are the saved derivative coefficients
!OUTPUT:
! grad_x, grad_y, grad_z = the x, y, and z components of the the gradient
!HISTORY:
! 18-may-2006, jmm, jimm@ssl.berkeley.edu
! 15-feb-2007, fixed bug which used cxarr for all deriv calls, jmm
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine grad_3d(a, x, y, z, nx, ny, nz, use_saved, 
     &     qsphere, grad_x, grad_y, grad_z, cxarr, cyarr, czarr)
      implicit none
      integer nx, ny, nz, use_saved, qsphere
      integer :: j, k, l
      real*8 :: a(nx, ny, nz), grad_x(nx, ny, nz), grad_y(nx, ny, nz),
     &     grad_z(nx, ny, nz)
      real*8 :: cxarr(nx, 3), cyarr(ny, 3), czarr(nz, 3)
      real*8 :: x(nx), y(ny), z(nz)
      real*8,allocatable::sinth(:)
      call deriv_dx(x, a, nx, ny, nz, use_saved, grad_x, cxarr)
      call deriv_dy(y, a, nx, ny, nz, use_saved, grad_y, cyarr)
      call deriv_dz(z, a, nx, ny, nz, use_saved, grad_z, czarr)
      if(qsphere.eq.1) then
         allocate(sinth(ny))
         sinth = sin(y)
         forall (j=1:nx,k=1:ny,l=1:nz)
            grad_y(j,k,l) = grad_y(j,k,l)/x(j)
            grad_z(j,k,l) = grad_z(j,k,l)/(x(j)*sinth(k))
         end forall
         deallocate(sinth)
      end if
      return
      end
!+
!NAME:
! obj_funct_fff
!PURPOSE:
! Calculates the objective function for the FFF code,
! See Wheatland, etal, 2000, eqn(5) Obj funct is
! the integral of b2*Omega2 over the whole volume.
!CALLING SEQUENCE:
! L = obj_funct_fff(dx, dy, dz, x, y, z, nx, ny, nz,
!                   b2, omega2, qsphere, use_saved)
!INPUT:
! dx, dy, dz, the differentials of the spatial grid
! x, y, z = the spatial grid, held in the bfields common block
! b2 = B2 
! omega2 =omega ^2 ; omega = (1/B2)*((curlB X B) - (div_B)B)
! qsphere = set to 1 for spherical coordinates
! use_saved =  set to 1, use saved grids
!OUTPUT:
! L = the integral of B2*Omega2 over the volume
!HISTORY:
! 19-May-2006, jmm
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      Function obj_funct_fff(x, y, z, nx, ny, nz, b2, 
     &     omega2, qsphere, use_saved, coeffx, coeffy, coeffz)
      implicit none

      integer :: nx, ny, nz, use_saved, qsphere
      integer :: j, k, l
      real*8 :: x(nx), y(ny), z(nz)
      real*8 :: omega2(nx, ny, nz), b2(nx, ny, nz)
      real*8 :: coeffx(nx), coeffy(ny), coeffz(nz)
      real*8 :: obj_funct_fff
      real*8,allocatable:: dx(:), dy(:), dz(:)
      if(use_saved.Eq.0) then
         allocate (dx(nx), dy(ny), dz(nz))
         do 5 j=1,nx-1
            dx(j)=x(j+1)-x(j)
 5       continue
         do 10 k=1,ny-1
            dy(k)=y(k+1)-y(k)
 10      continue
         do 15 l=1,nz-1
            dz(l)=z(l+1)-z(l)
 15      continue
!Arrays for integration
         coeffx(1)=0.5*dx(1)
         coeffx(nx)=0.5*dx(nx-1)
         coeffx(2:nx-1) = 0.5*(dx(1:nx-2)+dx(2:nx-1))
         coeffy(1)=0.5*dy(1)
         coeffy(ny)=0.5*dy(ny-1)
         coeffy(2:ny-1) =  0.5*(dy(1:ny-2)+dy(2:ny-1))
         coeffz(1)=0.5*dz(1)
         coeffz(nz)=0.5*dz(nz-1)
         coeffz(2:nz-1) =  0.5*(dz(1:nz-2)+dz(2:nz-1))
         If(qsphere.eq.1) Then
            coeffx = (x**2)*coeffx
            coeffy = sin(y)*coeffy
         Endif
         deallocate(dx, dy, dz)
      end if
      obj_funct_fff=0.0
      do 40 j=1,nx
         do 30 k=1, ny
            do 20 l=1,nz
               obj_funct_fff=obj_funct_fff+b2(j,k,l)*omega2(j,k,l)*
     &              coeffx(j)*coeffy(k)*coeffz(l)
 20         continue
 30      continue
 40   continue
      return
      end
!+
!NAME:
! calc_lin_fff
!PURPOSE:
! Calculation to determine linear FFF from the field on the lower
! boundary. Note that this only works for cartesian coordinates
!CALLING SEQUENCE:
! call calc_lin_fff(bx, by, bz, x, y, z, 
!     &     rsize, tsize, vsize, qsphere, alpha)
!INPUT:
! bx, by, bz = the magnetic field
! x, y, z = the spatial grid
!OUTPUT:
! bx, by, bz = the poential field solution, except on the lower
!              boundary, which is unchanged
!HISTORY:
! 1-jun-2006, jmm, jimm@ssl.berkeley.edu
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine Calc_lin_fff(bx, by, bz, x, y, z, 
     &     nx, ny, nz, qsphere, alpha)
      implicit none
      integer :: nx, ny, nz, qsphere
      integer :: j, k, l, jp, kp
      real*8 :: alpha
      real*8 :: x(nx), y(ny), z(nz)
      real*8 :: bx(nx, ny, nz), by(nx, ny, nz), bz(nx, ny, nz)
      real*8, allocatable:: ddx(:), ddy(:)
      real*8, allocatable:: dx(:), dy(:)
      real*8, allocatable:: cos_az(:), sin_az(:)
      real*8 r, bigr, cos_ar, sin_ar, xsum, ysum, zsum
      real*8 g, dgdz, gx, gy, gz, twopie, zl
  
      If(alpha.Ne.0.0) Then
         write(6,*) 'Linear FFF implemented, alpha = ', alpha
      Else 
         write(6,*) 'Potential Field implemented'
      End if
!Arrays for integration
      twopie = 2.0*3.1415927
      allocate (dx(nx), dy(ny))
      allocate (ddx(nx), ddy(ny))
      allocate (cos_az(nz), sin_az(nz))
      do 5 j=1,nx-1
         dx(j)=x(j+1)-x(j)
 5    continue
      do 10 k=1,ny-1
         dy(k)=y(k+1)-y(k)
 10   continue
      ddx(1)=0.5*dx(1)
      ddx(nx)=0.5*dx(nx-1)
      ddx(2:nx-1) = 0.5*(dx(1:nx-2)+dx(2:nx-1))
      ddy(1)=0.5*dy(1)
      ddy(ny)=0.5*dy(ny-1)
      ddy(2:ny-1) =  0.5*(dy(1:ny-2)+dy(2:ny-1))
!get cos and sin arrays
      cos_az = cos(alpha*z)
      sin_az = sin(alpha*z)
!5 loops, for each j,k,l, integrate over x and y
      do 60 l=2,nz
         zl = z(l)-z(1)
         do 50 j=1,nx
            do 40 k=1,ny
               xsum=0.0
               ysum=0.0
               zsum=0.0
               do 30 jp=1,nx 
                  do 20 kp=1,ny
                     if(jp.ne.j.or.kp.ne.k) then
                        bigR = sqrt((x(j)-x(jp))**2+(y(k)-y(kp))**2)
                        r = sqrt(bigR**2+zl**2)
                        cos_ar = cos(alpha*r)
                        sin_ar = sin(alpha*r)
                        g = zl*cos_ar/(bigR*r)-cos_az(l)/bigR
                        dgdz = cos_ar*(1.0/(bigr*r)-
     &                       zl**2/(bigr*r**3))-
     &                       alpha*zl**2*sin_ar/(bigr*r**2)+
     &                       alpha*sin_az(l)/bigr
                        gx = bz(jp,kp,1)*((x(j)-x(jp))*dgdz/bigr+
     &                       alpha*g*(y(k)-y(kp))/bigr)
                        gy = bz(jp,kp,1)*((y(k)-y(kp))*dgdz/bigr-
     &                       alpha*g*(x(j)-x(jp))/bigr)
                        gz = bz(jp,kp,1)*(zl*cos_ar/r**3+
     &                       alpha*zl*sin_ar/r**2)
!                        if(j.eq.10.and.k.eq.4.and.l.eq.2) then
!                           if(jp.eq.1.and.kp.eq.1) then
!                              write(6,*) x(j), y(k)
!                              write(6,*) x(jp), y(kp)
!                              write(6,*) bigr, r, cos_ar,sin_ar
!                              write(6,*) g, dgdz
!                              write(6,*) bz(jp,kp,1)
!                              write(6,*) gx, gy, gz
!                              write(6,*) ddx(jp), ddy(kp)
!                           endif
!                        endif
                        xsum = xsum+gx*ddx(jp)*ddy(kp)
                        ysum = ysum+gy*ddx(jp)*ddy(kp)
                        zsum = zsum+gz*ddx(jp)*ddy(kp)
                     end if
 20               continue
 30            continue
               bx(j,k,l)=xsum/twopie
               by(j,k,l)=ysum/twopie
               bz(j,k,l)=zsum/twopie
!               if(j.eq.10.and.k.eq.4.and.l.eq.2) then
!                  write(6,*) bx(j,k,l),by(j,k,l),bz(j,k,l)
!                  stop
!               endif
 40         continue
 50      continue
 60   continue

      deallocate (dx, dy)
      deallocate (ddx, ddy)
      deallocate (cos_az, sin_az)

      Return
      End
