subroutine load_quads (u, frac)

use cesrv_struct
use cesrv_interface
use q_tune_mod

implicit none

type (universe_struct), target :: u

real(rp) fh_fix, fv_fix, part, ky_fix_tot, kx_fix_tot, cu_tune_fix, frac
real(rp) dkx_dqa, dkx_dqb, dky_dqa, dky_dqb, f_khz
real(rp), save :: old_frac

integer, parameter :: MAX_CHANGE=4095
integer i, jj, ii, ix, save_set, ios, isteps, istat
integer cu0_csr_quad_cur(98), load_csr_quad_cur(98), target_csr_quad_cur(98)
integer cu0_csr_qadd_cur(98), load_csr_qadd_cur(98), target_csr_qadd_cur(98)
integer cu0_csr_sqewquad(98), load_csr_sqewquad(98), target_csr_sqewquad(98)
integer delt_csr_quad_cur(98), delt_csr_qadd_cur(98), delt_csr_sqewquad(98)
integer do_step, vnumbr
integer largest_change, substeps, n_quad
real(rp) subpart

character(80) line

logical limited, model_is_design, quad_model_is_design

! If the model = design then we are continuing form the last load.
! If so then transfer old -> model.
! The exception is If we are trying to load after reading in a data set.
! In this case we are trying to restore the saved
!   values and logic%use_old_in_loading = F

model_is_design = .true.
quad_model_is_design = .true.

do i = lbound(u%quad_k1%v, 1), ubound(u%quad_k1%v, 1)
  if (u%quad_k1%v(i)%good_var .and. u%quad_k1%v(i)%model /= u%quad_k1%v(i)%design) then
    model_is_design =  .false.
    quad_model_is_design =  .false.
    exit
  endif
enddo

do i = lbound(u%skew_quad_k1%v, 1), ubound(u%skew_quad_k1%v, 1)
  if (u%skew_quad_k1%v(i)%good_var .and. u%skew_quad_k1%v(i)%model /= u%skew_quad_k1%v(i)%design) then
    model_is_design =  .false.
    exit
  endif
enddo

if (model_is_design .and. logic%use_old_in_loading) then
  print *, 'Note: Looks like you are doing a second load...'
  print *, '      Reinstating model'
  call do_var_transfer (u%quad_k1, 1.0_rp, 'OLD', u, .false.)
  call do_var_transfer (u%skew_quad_k1, 1.0_rp, 'OLD', u, .false.)

! q_tune the model quads so the corrections will not change the tune?

else
  old_frac = 0
  n_quad = 0
  do i = lbound(u%quad_k1%v, 1), ubound(u%quad_k1%v, 1)
    if (.not. q_tune_use_this_quad(i, u, .true.)) cycle
    n_quad = n_quad + 1
  enddo

  ! Qtune?
  1000 continue
  if (.not. quad_model_is_design .and. (abs(u%tune%x%d(1)%meas - u%tune%x%d(1)%design) > 0.01 .or. &
                                        abs(u%tune%y%d(1)%meas - u%tune%y%d(1)%design) > 0.01_rp)) then

    if (n_quad == 0) then
      print *, 'NOTE: Model tunes do not match Design tunes so tunes will walk when the correction is loaded.'
      print *, '  Since there are no quadrupoles in optimization list there is nothing to do about this.'
    else
      print *, '[Note: q_tuneing here only uses the quads in the optimization list.]'
      print *, 'Change Model tunes to match Design so tunes will not'
      call get_input_string ('walk when you put in the corrections? <CR = Yes>', line)
      call string_trim (line, line, ix)
      call str_upcase(line, line)
      if (ix == 0 .or. index('YES', line(1:ix)) == 1)  then
        call q_tune (u%tune%x%d(1)%design, u%tune%y%d(1)%design, u, .true.)
      elseif (index('NO', line(1:ix)) /= 1) then
        print *, 'I DO NOT UNDERSTAND THIS. TRY AGAIN...'
        goto 1000
      endif
    endif

  endif

endif

! number of steps

isteps = 10
istat = in4get1 ('Number of steps to take? @', isteps)
if (istat /= 1) return

! check for quads out of range

call limit_calc (limited)    ! find cu values
if (limited) then
  print *, 'SINCE THE QUAD STRENGTH HAS SHIFTED BY THE LIMIT CALC, NOTHING WILL BE DONE'
  return
endif

! check the current save set

print *, 'QUADRUPOLE CORRECTION LOADING:'
print *

call vxgetn ('CSR SAVRECRD', 1, 1, save_set)  ! get save set number
if (save_set /= u%phase%csr_set) then
  print *, 'WARNING: YOU ARE NOT STARTING FROM THE SAVE SET IN WHICH THE'
  print *, '         DATA WAS TAKEN  [This is not necessarily an error. '
  print *, '         The end result is always the same irregardless.]'
  print *, '         Data save set:   ', u%phase%csr_set
  print *, '         Current save set:', save_set
  if (.not. logic%gui) then
  call get_input_string ('Continue? <CR = No>:', line)
    if (line(1:1) /= 'Y' .and. line(1:1) /= 'y') return
  endif
endif

! get current settings

ix = vnumbr('CSR QADD CUR')
call vxgetn ('CSR QUAD CUR', 1, 98, cu0_csr_quad_cur)
call vxgetn ('CSR QADD CUR', 1, ix, cu0_csr_qadd_cur)
call vxgetn ('CSR SQEWQUAD', 1, 98, cu0_csr_sqewquad)

! get target values

call cu_target_calc (u%var, u, frac+old_frac)

do i = 1, size(u%db%csr_quad_cur, 1)
  ix = u%db%csr_quad_cur(i)%ix_cesrv
  if (ix < 0) cycle
  target_csr_quad_cur(i) = u%quad_k1%v(ix)%cu_target
enddo

do i = 1, size(u%db%csr_qadd_cur, 1)
  ix = u%db%csr_qadd_cur(i)%ix_cesrv
  if (ix < 0) cycle
  target_csr_qadd_cur(i) = u%quad_k1%v(ix)%cu_target
enddo

do i = 1, size(u%db%csr_sqewquad, 1)
  ix = u%db%csr_sqewquad(i)%ix_cesrv
  if (ix < 0) cycle
  target_csr_sqewquad(i) = u%skew_quad_k1%v(ix)%cu_target
enddo

! get number of steps

if (isteps < 1) then
  print *, 'ERROR: NUMBER OF STEPS MUST BE POSITIVE.'
  print *, '       NOTHING LOADED'
  return
else
  print *, 'Loading steps:', isteps
endif

! put in changes

ky_fix_tot = 0. !total fractional change in vf quads
kx_fix_tot = 0. !total fractional change in hf quads
call q_tune_coef_calc (u, dkx_dqa, dkx_dqb, dky_dqa, dky_dqb, .true.)
f_khz = twopi / 390

do ii = 1, isteps

  part = float(ii) / isteps

  if (logic%gui) then    ! Get instructions from the gui
    call check_load_gui(do_step, fh_fix, fv_fix)
    if (do_step /= 1) exit    ! Quit loading
    print *,' fh_fix = ',fh_fix,' fv_fix = ',fv_fix
  else
    2000 continue
    call get_input_string (&
        'Enter change in tune (kHz) [fh, fv] (or "QUIT") <CR = 0, 0>:', line)
    call string_trim (line, line, ix)
    call str_upcase (line, line)
    if (ix == 0) then
      fh_fix = 0
      fv_fix = 0
    elseif (index('QUIT', line(:ix)) == 1) then
      frac = (ii - 1) * frac / isteps
      exit
    else
      read (line, *, iostat = ios) fh_fix, fv_fix
      if (ios /= 0) then
        print *
        print *, 'ERROR: CANNOT READ TUNE CHANGES. TRY AGAIN...'
        print *
        goto 2000
      endif
    endif
  endif

  kx_fix_tot = kx_fix_tot + f_khz * (dkx_dqa * fh_fix + dkx_dqb * fv_fix)
  ky_fix_tot = ky_fix_tot + f_khz * (dky_dqa * fh_fix + dky_dqb * fv_fix)

  ! Calculate quad changes

  do jj = 1, 98

    if (.not. q_tune_use_this_quad (jj, u, .true.)) then
      cu_tune_fix = 0
    elseif (u%quad_k1%v(jj)%design < 0) then    ! if vert focusing
      cu_tune_fix = ky_fix_tot / u%quad_k1%v(jj)%dvar_dcu
    else
      cu_tune_fix = kx_fix_tot / u%quad_k1%v(jj)%dvar_dcu
    endif

    delt_csr_quad_cur(jj) = nint(part * &
         (target_csr_quad_cur(jj) - cu0_csr_quad_cur(jj)) + cu_tune_fix)

    delt_csr_qadd_cur(jj) = nint(part * &
         (target_csr_qadd_cur(jj) - cu0_csr_qadd_cur(jj)))

    delt_csr_sqewquad(jj) = nint(part * &
         (target_csr_sqewquad(jj) - cu0_csr_sqewquad(jj)))

  enddo

  ! find the total changes to be loaded, and how many steps required
  ! to keep the change for each step < 4095
  largest_change = 0
  largest_change = max(maxval(abs(delt_csr_quad_cur)), largest_change)
  largest_change = max(maxval(abs(delt_csr_qadd_cur)), largest_change)
  largest_change = max(maxval(abs(delt_csr_sqewquad)), largest_change)

  substeps = ceiling(float(largest_change) / MAX_CHANGE)

  ! make sure scaler values match commands
  call set_sca_to_cmd

  ! put in changes

  call vcmode ('shortd')
  do jj=1,substeps
    subpart = float(jj) / substeps
    load_csr_quad_cur = nint(cu0_csr_quad_cur + subpart * delt_csr_quad_cur)
    load_csr_qadd_cur = nint(cu0_csr_qadd_cur + subpart * delt_csr_qadd_cur)
    load_csr_sqewquad = nint(cu0_csr_sqewquad + subpart * delt_csr_sqewquad)

    if (logic%debug) then
      print *, 'Fake Loading: QUADS and SKEW_QUADS'
      print '(a, (/, 5x, 10i6))', 'CSR QUAD CUR', load_csr_quad_cur(1:98)
      print '(a, (/, 5x, 10i6))', 'CSR QADD CUR', load_csr_qadd_cur(1:n_qadd_maxx)
      print '(a, (/, 5x, 10i6))', 'CSR SKEWQUAD', load_csr_sqewquad(1:98)
    else
      call vxputn ('CSR QUAD CUR', 1, 98,          load_csr_quad_cur)
      ix = vnumbr('CSR QADD CUR')
      call vxputn ('CSR QADD CUR', 1, ix, load_csr_qadd_cur)
      call vxputn ('CSR SQEWQUAD', 1, 98,          load_csr_sqewquad)
      call vclogo
    endif

  enddo
  call vcfree
  call vcmode ('')   ! Put back to immediate mode.

  print '(a, i3)', ' Quad k1 changes loaded for step: ', ii

enddo

! zero model

do i = lbound(u%quad_k1%v, 1), ubound(u%quad_k1%v, 1)
  u%quad_k1%v(i)%old = u%quad_k1%v(i)%model
  u%quad_k1%v(i)%model = u%quad_k1%v(i)%design
  call var_bookkeeper (u%quad_k1%v(i), u%ring, u%orb)
enddo

do i = lbound(u%skew_quad_k1%v, 1), ubound(u%skew_quad_k1%v, 1)
  u%skew_quad_k1%v(i)%old = u%skew_quad_k1%v(i)%model
  u%skew_quad_k1%v(i)%model = u%skew_quad_k1%v(i)%design
  call var_bookkeeper (u%skew_quad_k1%v(i), u%ring, u%orb)
enddo

print *, "Note: MODEL Quad K1's and Skew_quad K1's have been reset to DESIGN"
print *, "      [MODEL values are saved in OLD.]"

logic%use_old_in_loading = .true.
print '(a, f8.2)', 'Quadrupoles loaded. Total fraction loaded:', frac+old_frac
old_frac = frac + old_frac

end subroutine