!+
!NAME:
! evolve_fff
!PURPOSE:
! Evolves the magnetc field for 1 time step
!CALLING SEQUENCE:
! call evolve_fff(bx, by, bz, wf, x, y, z, 
!         nx, nz, ny, dt, qsphere, imax 
!         abs_frac_diff, local_l, conv_flag)
!INPUT:
! bx, by, bz = the magnetic field (also output)
! wf = weighting function
! x, y, z, = the spatial grid
! nx, nz, ny = the grid sizes
! dt = step size in "Time" t
! qsphere = set to 1 for spherical coordinates
! imax = the max number of iterations
! abs_frac_diff = the fractional change in the l_value needed
!                 for convergence (the change in L needs to be 
!                 smaller than this).
!OUTPUT:
! local_l = final value of l_value
! conv_flag = set to 1 if converged..
!HISTORY:
! 22-may-2006, jmm, jimm@ssl.berkeley.edu
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
!-
      subroutine evolve_fff(bx, by, bz, wf, x, y, z, 
     &     nx, ny, nz, qsphere, imax, 
     &     dt, abs_frac_diff, local_l, conv_flag)
      implicit none
 
      integer :: nx, ny, nz, use_saved, i, imax
      integer :: nx1, ny1, nz1, qsphere, conv_flag
      
      real*8 :: dt, l_value, local_l, abs_frac_diff, delta_l, dtmin
      real*8 :: x(nx), y(ny), z(nz)
      real*8 :: bx(nx, ny, nz), by(nx, ny, nz), bz(nx, ny, nz)
      real*8 :: wf(nx, ny, nz)
      real*8, external :: obj_funct_fff

      real*8, allocatable:: coeffx(:), coeffy(:), coeffz(:)
      real*8, allocatable:: cxarr(:,:), cyarr(:,:), czarr(:,:)

      real*8, allocatable:: jx(:,:,:), jy(:,:,:), jz(:,:,:)
      real*8, allocatable:: omega_x(:,:,:), omega_y(:,:,:), 
     &     omega_z(:,:,:)
      real*8, allocatable:: div_b(:,:,:), omega2(:,:,:), b2(:,:,:)
      real*8, allocatable:: f_x(:,:,:), f_y(:,:,:), f_z(:,:,:)

      allocate(coeffx(nx),coeffy(ny),coeffz(nz))
      allocate(cxarr(nx,3),cyarr(ny,3),czarr(nz,3))
      allocate(jx(nx,ny,nz),jy(nx,ny,nz),jz(nx,ny,nz))
      allocate(omega_x(nx,ny,nz),omega_y(nx,ny,nz),omega_z(nx,ny,nz))
      allocate(div_b(nx,ny,nz),omega2(nx,ny,nz),b2(nx,ny,nz))
      allocate(f_x(nx,ny,nz),f_y(nx,ny,nz),f_z(nx,ny,nz))      

!need some parameters
      nx1=nx-1
      nz1=nz-1
      ny1=ny-1
      i=0
      conv_flag=0
      l_value=1.0e20            !high number to start...
      delta_l=1.0e20
      dtmin=1.0e-10
      open(unit=18,file='lhist.dat',status='replace')
      do while (i.lt.imax.and.conv_flag.eq.0.and.dt.gt.dtmin)
         i=i+1
         if(i.eq.1) then        !on the first iteration, i need this
            use_saved=0
            b2 = bx**2+by**2+bz**2
            call curl_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &           use_saved, qsphere, jx, jy, jz, cxarr, cyarr, czarr)
            call div_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &           use_saved, qsphere, div_b, cxarr, cyarr, czarr)
! 20      format(5f16.7)
!         open(unit=8,file='field_test_b.dat',status='replace')
!         write(8,*) nx, ny, nz
!         write(8,20) bx
!         write(8,20) by
!         write(8,20) bz
!         write(8,20) x
!         write(8,20) y
!         write(8,20) z
!         close(unit=8)
!         open(unit=8,file='field_test_j.dat',status='replace')
!         write(8,*) nx, ny, nz
!         write(8,20) jx
!         write(8,20) jy
!         write(8,20) jz
!         write(8,20) x
!         write(8,20) y
!         write(8,20) z
!         close(unit=8)
!         open(unit=8,file='field_test_divb.dat',status='replace')
!         write(8,*) nx, ny, nz
!         write(8,20) div_b
!         close(unit=8)
!        stop
            omega_x = jy*bz-jz*by
            omega_y = jz*bx-jx*bz
            omega_z = jx*by-jy*bx
            where(b2.ne.0) 
               omega_x = (omega_x - div_b*bx)/b2
               omega_y = (omega_y - div_b*by)/b2
               omega_z = (omega_z - div_b*bz)/b2
            elsewhere
               omega_x = 0.0
               omega_y = 0.0
               omega_z = 0.0
            end where
            omega2 = wf*(omega_x**2+omega_y**2+omega_z**2)
            omega_x = omega_x*wf
            omega_y = omega_y*wf
            omega_z = omega_z*wf
            l_value = obj_funct_fff(x, y, z, nx, ny, nz, 
     &           b2, omega2, qsphere, use_saved,
     &           coeffx, coeffy, coeffz)
            use_saved=1
         end if                 !use omega and B that you have
!         if(mod(i,10).eq.1) then
            write(18,*) i, dt, l_value
            write(6,*) i, dt, l_value, delta_l
!         endif
!now get the force, 
!maybe we can save memory, subtract omega cross j
         f_x = -(omega_y*jz-omega_z*jy)
         f_y = -(omega_z*jx-omega_x*jz)
         f_z = -(omega_x*jy-omega_y*jx)
!Now you can use b2 and jx,jy,jz for temp storage
!subtract grad omega dot B, 
!Add omega times div_b !Add omega2 times B
         b2 = omega_x*bx+omega_y*by+omega_z*bz
         call grad_3d(b2, x, y, z, nx, ny, nz, use_saved,
     &        qsphere, jx, jy, jz, cxarr, cyarr, czarr)
         f_x = f_x-jx+omega_x*div_b+omega2*bx
         f_y = f_y-jy+omega_y*div_b+omega2*by
         f_z = f_z-jz+omega_z*div_b+omega2*bz
! Add curl of omega cross b, b2 will be the x component,
! div_b will be the y component, omega2 the z component, 
!jx,jy,jz are the output
         b2 = omega_y*bz-omega_z*by
         div_b = omega_z*bx-omega_x*bz
         omega2 = omega_x*by-omega_y*bx
         call curl_3d(b2, div_b, omega2, x, y, z, 
     &        nx, ny, nz, use_saved, qsphere, 
     &        jx, jy, jz, cxarr, cyarr, czarr)
         f_x = f_x+jx
         f_y = f_y+jy
         f_z = f_z+jz
!new B
         bx(2:nx1,2:ny1,2:nz1) = bx(2:nx1,2:ny1,2:nz1)+
     &        dt*f_x(2:nx1,2:ny1,2:nz1)
         by(2:nx1,2:ny1,2:nz1) = by(2:nx1,2:ny1,2:nz1)+
     &        dt*f_y(2:nx1,2:ny1,2:nz1)
         bz(2:nx1,2:ny1,2:nz1) = bz(2:nx1,2:ny1,2:nz1)+
     &        dt*f_z(2:nx1,2:ny1,2:nz1)
!now new values for j, omega
         b2 = bx**2+by**2+bz**2
         call curl_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &        use_saved, qsphere, jx, jy, jz, cxarr, cyarr, czarr)
         call div_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &        use_saved, qsphere, div_b, cxarr, cyarr, czarr)
         omega_x = jy*bz-jz*by
         omega_y = jz*bx-jx*bz
         omega_z = jx*by-jy*bx
         where(b2.ne.0) 
            omega_x = (omega_x - div_b*bx)/b2
            omega_y = (omega_y - div_b*by)/b2
            omega_z = (omega_z - div_b*bz)/b2
         elsewhere
            omega_x = 0.0
            omega_y = 0.0
            omega_z = 0.0
         end where
         omega2 = wf*(omega_x**2+omega_y**2+omega_z**2)
         omega_x = omega_x*wf
         omega_y = omega_y*wf
         omega_z = omega_z*wf
         local_l = obj_funct_fff(x, y, z, nx, ny, nz, 
     &        b2, omega2, qsphere, use_saved,
     &        coeffx, coeffy, coeffz)
         delta_l = (l_value-local_l)/l_value
         if(delta_l.Gt.0) then !success
            dt = dt*1.01
!convergence?
            l_value = local_l
            if(abs(delta_l).lt.abs_frac_diff) then 
               conv_flag=1
               write(6,*) 'conv_flag = ', conv_flag
            endif
         else                   !reset B, j, and omega
            dt = dt/2.0         !reset dt
!reset b
            bx(2:nx1,2:ny1,2:nz1) = bx(2:nx1,2:ny1,2:nz1)-
     &           dt*f_x(2:nx1,2:ny1,2:nz1)
            by(2:nx1,2:ny1,2:nz1) = by(2:nx1,2:ny1,2:nz1)-
     &           dt*f_y(2:nx1,2:ny1,2:nz1)
            bz(2:nx1,2:ny1,2:nz1) = bz(2:nx1,2:ny1,2:nz1)-
     &           dt*f_x(2:nx1,2:ny1,2:nz1)
!now new values for j, omega, b2, div_b, omega2
            b2 = bx**2+by**2+bz**2
            call curl_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &           use_saved, qsphere, jx, jy, jz, cxarr, cyarr, czarr)
            call div_3d(bx, by, bz, x, y, z, nx, ny, nz, 
     &           use_saved, qsphere, div_b, cxarr, cyarr, czarr)
            omega_x = jy*bz-jz*by
            omega_y = jz*bx-jx*bz
            omega_z = jx*by-jy*bx
            where(b2.ne.0) 
               omega_x = (omega_x - div_b*bx)/b2
               omega_y = (omega_y - div_b*by)/b2
               omega_z = (omega_z - div_b*bz)/b2
            elsewhere
               omega_x = 0.0
               omega_y = 0.0
               omega_z = 0.0
            end where
            omega2 = wf*(omega_x**2+omega_y**2+omega_z**2)
            omega_x = omega_x*wf
            omega_y = omega_y*wf
            omega_z = omega_z*wf
         end if
      end do
      write(18,*) i, dt, l_value
      write(6,*) i, dt, l_value, delta_l
      close(unit=18)
      if(conv_flag.eq.1) then 
         write(6,*) 'Successful convergence'
      else
         write(6,*) 'Not converged'
      endif
      deallocate(coeffx,coeffy,coeffz)
      deallocate(cxarr,cyarr,czarr)
      deallocate(jx,jy,jz)
      deallocate(omega_x,omega_y,omega_z)
      deallocate(div_b,omega2,b2)

      return

      end
