!==============================================================!
! PROGRAM optimization_fff                                     !
!==============================================================!
! PURPOSE:                                                     !
!   Program to extrapolate linear force-free fields, based on  !
!   the optimization algorithm published in Apj 540, 1150.     !
! INPUT:                                                       !
!   3-d initial field and grids (e.g., as output by            !  
!   write_bfield_fff.pro) in file 'field0.dat'                 !
!                                                              !
!   Other parameters are read in from a file called            !
!   'fff_input_pars.dat', which contains tha maximum number of !
!   iterations (imax), the spherical coordinates logical flag  !
!   (qsphere), the convergence parameter and the initial value !
!   of dt.                                                     !
! OUTPUT:                                                      !
!   The evolved field in the file 'field.dat'                  !
!                                                              !
!   Diagnostic file 'lhist.dat'                                !
!==============================================================!
PROGRAM optimization_fff
  USE params
  USE vars
  USE io
  USE fff_operators
  USE fff_evolve
  IMPLICIT NONE
 
  ! Local variable declarations
  INTEGER(KIND=i4) :: fdata, nd, nbpw
  INTEGER(KIND=i8) :: d1
  CHARACTER(LEN=1) :: dtyp
  CHARACTER(LEN=10) :: lbl
  REAL(KIND=r8) :: bzmaxinv

  CALL openfile(ferror,'error.log')
  ! Read input parameters
  CALL read_input
  ! Read in data
#if SDF
  CALL sdf_read_f77('field0.sdf',0,lbl,dtyp,nbpw,nd,d1,nx)
  CALL sdf_read_f77('field0.sdf',1,lbl,dtyp,nbpw,nd,d1,ny)
  CALL sdf_read_f77('field0.sdf',2,lbl,dtyp,nbpw,nd,d1,nz)
  ALLOCATE (x(nx), y(ny), z(nz))
  ALLOCATE (bx(nx,ny,nz), by(nx,ny,nz), bz(nx,ny,nz))
  ALLOCATE (wf(nx,ny,nz))
  CALL read_sdf_data
#else
  CALL openfile(fdata,'field0.dat',IFORM='formatted',ISTATUS='old')
  READ(UNIT=fdata,FMT=*) nx, ny, nz
  CALL closefile(fdata)
  ALLOCATE (x(nx), y(ny), z(nz))
  ALLOCATE (bx(nx,ny,nz), by(nx,ny,nz), bz(nx,ny,nz))
  ALLOCATE (wf(nx,ny,nz))
  CALL read_data
#endif

  ! Normalize the field by max value of bz(:,:,1)
  IF(qsphere .EQV. .TRUE.) THEN
    bzmax = MAXVAL(bz(1, :, :))
  ELSE
    bzmax = MAXVAL(bz(:, :, 1))
  END IF
  bzmaxinv = 1.0_r8/bzmax
  bx = bzmaxinv*bx
  by = bzmaxinv*by
  bz = bzmaxinv*bz
  ! If slow is set to true, redo the calculation using the Greens function
  IF(slow .EQV. .TRUE.) THEN
    CALL calc_lin_fff
    ! unnormalize the field by max value of bz(:,:,1)
    bx = bzmax*bx
    by = bzmax*by
    bz = bzmax*bz
#if SDF
    CALL write_sdf_field('ifield.sdf')
#else
    CALL write_field('ifield.dat')
#endif
    bx = bzmaxinv*bx
    by = bzmaxinv*by
    bz = bzmaxinv*bz
  ENDIF
  ! Evolve the field
  CALL evolve_fff
  ! Unnormalize the field by max value of bz(:,:,1)
  bx = bzmax*bx
  by = bzmax*by
  bz = bzmax*bz
#if SDF
  CALL write_sdf_field('field.sdf')
#else
  CALL write_field('field.dat')
#endif
  DEALLOCATE(x,y,z,bx,by,bz,wf)
  CALL shutdown('optimization_fff: clean exit')

END PROGRAM optimization_fff
