!+
!NAME:
!   optimization_fff
!PURPOSE:
!   Program to extrapolate linear force-free fields, based on the
!   optimization algorithm published in Apj 540, 1150.
!CALLING SEQUENCE:
!   This is a main program, compile by:
!   ifort -static-libcxa -o optimization_fff.exe fff_derivs.f 
!               fff_operators.f evolve_fff.f optimization_fff.f
!INPUT:
!   The input must be a 3-d initial field and grids, as output by 
!   write_bfield_fff.pro in the file 'field0.dat'
!   Other parameters are read in from a file called 
!   'fff_input_pars.dat', which contains the maximum number of 
!   iterations (imax), the spherical coordinates identifier
!   (1 for spherical, 0 for cartesian), the convergence 
!   parameter, and the initial value of dt.
!OUTPUT:
!   The field in the file 'field.dat'
!   Also a file 'lhist.dat' is created
!HISTORY:
!   2006, jmm, jimm@ssl.berkeley.edu
! 2-mar-2007, implicit none, and assignment statements added by Dave Bercik
! 14-sep-2007, added read - write from SDF files
!-   
      program optimization_fff
      implicit none
 
      integer :: rsize, tsize, vsize
      integer :: qsphere, conv_flag, imax, slow, sdfu
      integer :: j, k, l, nbpw, nd
      integer(kind=8) :: drstsvs(3), d1
      real*8 :: abs_frac_diff, dt, bnorm, alpha, local_l
      real*8, allocatable :: x(:), y(:), z(:)
      real*8, allocatable :: bx(:,:,:), by(:,:,:), bz(:,:,:)
      real*8, allocatable :: wf(:,:,:)
      character(len=10) lbl
      character(len=1) dtyp 
!first read the input variables
      open(unit=12,file='fff_input_pars.dat',status='old')
      read(12,*) dt
      read(12,*) abs_frac_diff
      read(12,*) imax
      read(12,*) qsphere
      read(12,*) slow
      read(12,*) alpha
      read(12,*) sdfu
      close(unit=12)
      If(sdfu.eq.1) Then
         call sdf_read_f77('field0.sdf',0,lbl,dtyp,nbpw,nd,d1,rsize)
         call sdf_read_f77('field0.sdf',1,lbl,dtyp,nbpw,nd,d1,tsize)
         call sdf_read_f77('field0.sdf',2,lbl,dtyp,nbpw,nd,d1,vsize)
         allocate (x(rsize), y(tsize), z(vsize))
         allocate (bx(rsize,tsize,vsize), by(rsize,tsize,vsize), 
     &        bz(rsize,tsize,vsize))
         call sdf_read_f77('field0.sdf',3,lbl,dtyp,nbpw,nd,drstsvs,bx)
         call sdf_read_f77('field0.sdf',4,lbl,dtyp,nbpw,nd,drstsvs,by)
         call sdf_read_f77('field0.sdf',5,lbl,dtyp,nbpw,nd,drstsvs,bz)
         call sdf_read_f77('field0.sdf',6,lbl,dtyp,nbpw,nd,d1,x)
         call sdf_read_f77('field0.sdf',7,lbl,dtyp,nbpw,nd,d1,y)
         call sdf_read_f77('field0.sdf',8,lbl,dtyp,nbpw,nd,d1,z)
         call sdf_read_f77('weight0.sdf',0,lbl,dtyp,nbpw,nd,d1,rsize)
         call sdf_read_f77('weight0.sdf',1,lbl,dtyp,nbpw,nd,d1,tsize)
         call sdf_read_f77('weight0.sdf',2,lbl,dtyp,nbpw,nd,d1,vsize)
         allocate (wf(rsize,tsize,vsize))
         call sdf_read_f77('weight0.sdf',3,lbl,dtyp,nbpw,nd,drstsvs,wf)
      Else
         open(unit=9,file='weight0.dat',status='old')
         read(9,*) rsize, tsize, vsize
         allocate (wf(rsize,tsize,vsize))
         read(9,*) wf
         close(unit=9)
         open(unit=11,file='field0.dat',status='old')
         read(11, *) rsize, tsize, vsize
         allocate (x(rsize), y(tsize), z(vsize))
         allocate (bx(rsize,tsize,vsize), by(rsize,tsize,vsize), 
     &        bz(rsize,tsize,vsize))
 10      format(3i12)
 20      format(5f16.7)
         read(11,*) bx
         read(11,*) by
         read(11,*) bz
         read(11,*) x
         read(11,*) y
         read(11,*) z
         close(unit=11)
      End if
!normalize the field by max value of bz(*,*,1)
      If(qsphere.eq.1) Then
         bnorm = maxval(bx(1, 1:tsize, 1:vsize))
      else
         bnorm = maxval(bz(1:rsize, 1:tsize, 1))
      end if
      forall(j=1:rsize,k=1:tsize,l=1:vsize)
         bx(j,k,l)=bx(j,k,l)/bnorm
         by(j,k,l)=by(j,k,l)/bnorm
         bz(j,k,l)=bz(j,k,l)/bnorm
      end forall
!if slow is 1, redo the calculation using the Greens function
      if(slow.eq.1) then
         call calc_lin_fff(bx, by, bz, x, y, z, 
     &     rsize, tsize, vsize, qsphere, alpha)
!unnormalize the field by max value of bz(*,*,1)
         forall(j=1:rsize,k=1:tsize,l=1:vsize)
            bx(j,k,l)=bx(j,k,l)*bnorm
            by(j,k,l)=by(j,k,l)*bnorm
            bz(j,k,l)=bz(j,k,l)*bnorm
         end forall
         open(unit=14,file='ifield.dat',status='replace')
         write(14,10) rsize, tsize, vsize
         write(14,20) bx
         write(14,20) by
         write(14,20) bz
         write(14,20) x
         write(14,20) y
         write(14,20) z
         close(unit=14)
         forall(j=1:rsize,k=1:tsize,l=1:vsize)
            bx(j,k,l)=bx(j,k,l)/bnorm
            by(j,k,l)=by(j,k,l)/bnorm
            bz(j,k,l)=bz(j,k,l)/bnorm
         end forall
      endif
!evolve the field
      call evolve_fff(bx, by, bz, wf, x, y, z, 
     &     rsize, tsize, vsize, qsphere, imax, 
     &     dt, abs_frac_diff, local_l, conv_flag)

!unnormalize the field by max value of bz(*,*,1)
      forall(j=1:rsize,k=1:tsize,l=1:vsize)
         bx(j,k,l)=bx(j,k,l)*bnorm
         by(j,k,l)=by(j,k,l)*bnorm
         bz(j,k,l)=bz(j,k,l)*bnorm
      end forall
      If(sdfu.eq.1) Then
         d1=1
         call sdf_write_f77('field.sdf','rsize','i',4,1,d1,rsize)
         call sdf_write_f77('field.sdf','tsize','i',4,1,d1,tsize)
         call sdf_write_f77('field.sdf','vsize','i',4,1,d1,vsize)
         drstsvs(1)=int(rsize,8)
         drstsvs(2)=int(tsize,8)
         drstsvs(3)=int(vsize,8)
         call sdf_write_f77('field.sdf','bx','f',8,3,drstsvs,bx)
         call sdf_write_f77('field.sdf','by','f',8,3,drstsvs,by)
         call sdf_write_f77('field.sdf','bz','f',8,3,drstsvs,bz)
         call sdf_write_f77('field.sdf','x','f',8,1,drstsvs(1),x)
         call sdf_write_f77('field.sdf','y','f',8,1,drstsvs(2),y)
         call sdf_write_f77('field.sdf','z','f',8,1,drstsvs(3),z)
      Else
         open(unit=10,file='field.dat',status='replace')
         write(10,10) rsize, tsize, vsize
         write(10,20) bx
         write(10,20) by
         write(10,20) bz
         write(10,20) x
         write(10,20) y
         write(10,20) z
         close(unit=10)
      End if

      End
