subroutine find_change(line, lat)

  use bmad_struct
  use bmad_interface
  use bmadz_interface

  implicit none



  type(lat_struct) lat
  type(lat_struct), save :: lat_0, lat_in, lat_out
  type(ele_struct) ele
  type (ele_struct), pointer :: slave
  type(coord_struct), allocatable :: co(:)
  type (all_pointer_struct) a_ptr
  type (ele_pointer_struct), allocatable :: eles(:), slave_eles(:)
  type (control_struct), pointer :: ctl

  real(rp) value
  real(rp), allocatable :: dk1(:)
  real(rp) phi_x, phi_y, int_Q_x, int_Q_y, qx, qy, qpx, qpy, Qz
  real(rp) val
  real(rp) current
  real(rp) energy
  real(rp) strength
  real(rp) mat(6,6)

  character*16 attribute
  character  line*(*)
  character *1 ans/'?'/
  character*16 element1, element2
  character*12 strength_word

  integer ix, n, ix_ele, ix_attr, change
  integer delta$/0/, setequal$/1/, proportional$/2/
  integer ele_key
  integer ixx, is
  integer i
  integer j,k 
  integer tot
  integer attrib
  integer particle, i_train, j_car, n_trains_tot, n_cars
  integer kee1, kee2
  integer lun
  integer n_loc
  integer number_elements/0/

  logical ok
  logical rec_taylor/.true./
  logical type_taylor/.false./

  logical, save :: first = .true.
  logical match
  logical err_flag
  logical err

!

  ix = index(line, '!')
  if (ix /= 0) line = line(:ix-1)        ! strip off comments


  if(first)lat_0 = lat
  first=.false.

  call str_upcase(line, line)
  call string_trim(line, line, ix)
 

    if(line(1:2) .eq. 'HE' .or. line(1:1) .eq. '?')then
 print '(a81)',' Type element name (Q04W), <= or delta or proportional>,<attribute>, and <value>'
 print '(a)','      Example 1:   Q04W = k1 0.2  (Result : Q04W[k1]=0.2)'
 print '(a)', &
 '      Example 2:   Q04W delta k1 0.001  (Result : Q04W[k1]=Q04W[k1]+0.001)'
 print '(a)', &
 '      Example 3:   Q04W proportional k1 1.1  (Result : Q04W[k1] =Q04W[k1]*1.1)'
  print '(a)', &
 '      Example 4:   H_SEP_08W delta hkick 0.0001'
 print '(a)','      Example 5:   SEX_08W delta k2 0.01'
 print '(a)','      Example 6:   RAW_QTUNEING_5 delta command 0.0001'
 print *
 print '(a)',' Hit <carriage return> and your last command will be repeated'
 print '(a)',' type "SHOW GROUP" to get a list of defined groups'
 print '(a)',' type "SHOW SKEW" to get a list of SKEW QUADS'
 print '(a)',' type "SHOW QUAD, <MARKER>,<DIPOLE>" to get a list of  QUADS, <MARKERS>,<DIPOLES>,<SEXTUPOLES>'
 print '(a)',' type "SHOW ELE <name>" to get details of element'
 print '(a)',' type "SHOW LAT <ix> <number>" to list elements and locations near in the range ix+- number'
 print '(a)',' type "SHOW TWISS <ix>" to list elements and twiss near ix'
 print '(a)',' type "SHOW MATRIX <ix>" to print 4X4 matrix at end of element ix (ix=0 for full turn)'
 print '(a)',' type "Q_TUNE  qx  qy" to qtune. type "Q_TUNE qx qy MATCH" to qtune using matching element'
 print '(a)',' type "QP_TUNE  qpx  qpy s1 s2" to set chromaticity'
 print '(a)',' type "QZ_TUNE  Qz" to set synchrotron tune'
 print '(a)',' type "ALL H_SEPARATOR value" to multiply all horizontal separators by <value>'
 print '(a)',' type "ALL V_SEPARATOR value" to multiply all vertical separators by <value>'
 print '(a)',' type "ALL SEXTUPOLES value" to multiply all sextupoles by <value>'
 print '(a)',' type "BEGINNING" to set twiss parameters at start'
 print '(a)',' type "LRBBI"  <particle> <i_train> <j_car> <n_trains_tot> <n_cars> <current(mA)> '
 print '(a)',' type "CBAR_V_E ON (OFF)" (see bmadz/code/cbar_v_e.f90'
 print '(a)',' type "SEXT_MOMENTS" to compute sexutpole moments'
 print '(a)','GO   (Result : implement changes and go)'
      return


  elseif(line(1:3) == 'SHO')then
    call string_trim(line(ix+1:),line,ix)
    if(line(1:3) == 'GRO')then
      do n=1,lat%n_ele_max
        if(index(lat%ele(n)%name,'RAW') /= 0)print *,n,lat%ele(n)%name
      end do
    endif
    if(line(1:3) == 'SKE')then
      do n=1,lat%n_ele_max
        if(index(lat%ele(n)%name,'SK') /= 0)print '(i10,1x,a16,es12.4)',n,lat%ele(n)%name, lat%ele(n)%value(k1$)
      end do
    endif
    
    kee1=0
    kee2=0
    if(line(1:3) == 'QUA')kee1=quadrupole$
    if(line(1:3) == 'MAR')kee1=marker$
    if(line(1:3) == 'SEX')kee1=sextupole$ 
    if(line(1:3) == 'WIG')kee1=wiggler$
    if(line(1:2) == 'RF')kee1=rfcavity$
    if(line(1:3) == 'DIP')then
     kee1=sbend$
     kee2=rbend$
    endif
    if(kee1 /= 0 .or. kee2 /= 0)then
     j=0

     do n=1,lat%n_ele_track
        if(kee1 == lat%ele(n)%key .or. kee2 == lat%ele(n)%key)then
             if(kee1 == quadrupole$)then
               strength = lat%ele(n)%value(k1$)
               strength_word = '    k1'
             endif
             if(kee1 == sextupole$)then
               strength = lat%ele(n)%value(k2$)
               strength_word = '    k2'
             endif
             if(line(1:3) == 'DIP')then
               strength = lat%ele(n)%value(b_field$)
               strength_word = '    B[T]'
             endif
             if(kee1== wiggler$)then
               strength = lat%ele(n)%value(b_max$)
               strength_word = '   B_max[T]'
             endif
             if(index(line,'&') == 0)then        
               if(j == 0)print '(a12,1x,a16,9a12)',' Element','Name            ','z',&
                                               'theta','beta x','beta y','eta x','alpha x','alpha y', &
                                                strength_word,'middle'
               print '(i12,1x,a16,9f12.4)',n,lat%ele(n)%name, lat%ele(n)%s, lat%ele(n)%floor%theta, &
                              lat%ele(n)%a%beta, lat%ele(n)%b%beta, lat%ele(n)%x%eta , &
                              lat%ele(n)%a%alpha, lat%ele(n)%b%alpha, strength, lat%ele(n)%s-lat%ele(n)%value(l$)/2
               j=j+1
               if(20*(j/20) == j .and. j /= 0)then
                 print '(a,$)',' Continue (C) or Exit (E) ?'
                 read(5,*)ans
                 call str_upcase(ans,ans)
                 if(ans(1:1) == 'E')exit
               endif

             else ! index(line,'&') /= 0        
                if(j == 0)then
                         lun = lunget()
                         open(unit = lun,file = line(1:3)//".dat")     
                         print '(a)',' write file '//line(1:3)//'.dat'
                         write(lun, '(a12,1x,a16,8a12)')' Element','Name            ','z',&
                                               'theta','beta x','beta y','eta x','alpha x','alpha y', &
                                                strength_word
                  endif
                    write(lun, '(i12,1x,a16,8f12.4)')n,lat%ele(n)%name, lat%ele(n)%s, lat%ele(n)%floor%theta, &
                         lat%ele(n)%a%beta, lat%ele(n)%b%beta, lat%ele(n)%x%eta , &
                         lat%ele(n)%a%alpha, lat%ele(n)%b%alpha, strength
                         j = j+1
                    if(n >= lat%n_ele_track)close(unit=lun)
                
            endif
        endif

      end do
    endif
    if(line(1:3) == 'ELE')then
      call string_trim(line(ix+1:),line,ix)
      call element_locator (line(1:ix), lat, ix_ele)
      print *
      if(ix_ele < 0)then
        print '(1x,a19,1x,a16)', ' No element named: ',line
       else
        call string_trim(line(ix+1:), line, ix)
        if(index(line,'TAYLOR')/= 0)type_taylor=.true.
        call type_ele(lat%ele(ix_ele), .false., 6, type_taylor, radians$)
        call type_v15(lat%ele(ix_ele))
       endif
    endif
    if(line(1:4) == 'LAT')then
      print *,line
      call string_trim(line(ix+1:),line,ix)
      print *,line
      if(line(1:2) == '  ')then
         print *,' SHO LAT <reference_element> <number_elements before and after>'
         return
      endif
      read(line(1:ix),*)ix_ele
      call string_trim(line(ix+1:),line,ix)
      number_elements = min(lat%n_ele_track,10)
      if(ix > 0)read(line(1:ix),*)number_elements
      print *
      print '(1x, a4,1x,a16,1x,9a14)',' ind',' element name   ','       s        ', &
         '      x         ','      y         ','      z        ','   theta       ','k1','k2','G','length'
      do i = ix_ele - min(number_elements, lat%n_ele_track), ix_ele + min(number_elements,lat%n_ele_track)
      j = i
      if(i < 0)j = lat%n_ele_track + i
      if(i > lat%n_ele_track)j = i - lat%n_ele_track
      print '(1x,i4,1x,a16,1x,9f14.6)', j, lat%ele(j)%name, lat%ele(j)%s, &
         lat%ele(j)%floor%r(1), lat%ele(j)%floor%r(2), &
         lat%ele(j)%floor%r(3), lat%ele(j)%floor%theta,lat%ele(j)%value(k1$),lat%ele(j)%value(k2$),lat%ele(j)%value(g$),lat%ele(j)%value(l$)
      end do
    endif
    if(line(1:4) == 'TWIS')then
      call string_trim(line(ix+1:),line,ix)
      read(line(1:ix),*)ix_ele
      print *
      print '(1x, a4,1x,a16,1x,6a12)','ind','element name','s', &
         'beta x','phase x','beta y',',phase y','eta x'
      do i = ix_ele - 10, ix_ele + 10
      j = i
      if(i < 0)j = lat%n_ele_track + i
      if(i > lat%n_ele_track)j = i - lat%n_ele_track
      print '(1x,i4,1x,a16,1x,6f12.4)', j, lat%ele(j)%name, lat%ele(j)%s, &
         lat%ele(j)%a%beta, lat%ele(j)%a%phi, &
         lat%ele(j)%b%beta, lat%ele(j)%b%phi, lat%ele(j)%x%eta
      end do
    endif
    if(line(1:4) == 'MATR')then
      call string_trim(line(ix+1:),line,ix)
      read(line(1:ix),*)ix_ele
      if(ix_ele <= 0 )ix_ele = lat%n_ele_track
      print *
      call mat_make_unit(mat)
      do i=1,ix_ele
       mat(1:6,1:6)= matmul(lat%ele(i)%mat6(1:6,1:6), mat(1:6,1:6))
      end do

      print '(a,a,a,i0)',' Transfer matrix from start to ', lat%ele(ix_ele)%name,' element ', ix_ele           
!      print '(1x, a4,1x,a16,1x,6a12)','ix_ele','element name','s'
     do i=1,6
        print '(6e12.4)',mat(i,1:6)
     end do
    endif

   elseif(line(1:3) == 'Q_T')then
     call string_trim(line(ix+1:), line, ix)
     read(line(1:ix),*)qx
     call string_trim(line(ix+1:), line, ix)
     if(ix /= 0)read(line(1:ix),*)qy
     call string_trim(line(ix+1:), line, ix)
     match = .false.
     if(ix /= 0)match = .true.

     call twiss_at_start(lat)
     call twiss_propagate_all(lat)
     int_Q_x = int(lat%ele(lat%n_ele_track)%a%phi / twopi)
     int_Q_y = int(lat%ele(lat%n_ele_track)%b%phi / twopi)
     phi_x = (int_Q_x + qx)*twopi
     phi_y = (int_Q_y + qy)*twopi


     allocate(co(0:lat%n_ele_max))
     allocate(dk1(lat%n_ele_max))
     call choose_quads(lat, dk1)
     forall(i=0:lat%n_ele_max)co(i)%vec(1:6)=0.
     call custom_set_tune (phi_x, phi_y, dk1, lat, co, ok, match= match)
     deallocate(co)
     deallocate(dk1)
     if (.not. ok)print *,' Qtune failed'

   elseif(line(1:3) == 'QP_')then
     call string_trim(line(ix+1:),line,ix)
     read(line(1:ix),*)qpx
     call string_trim(line(ix+1:),line,ix)
     read(line(1:ix),*)qpy

     element1 = 'RAW_XQUNEING_1'
     element2 = 'RAW_XQUNEING_2'
     call string_trim(line(ix+1:),line, ix)
     if(ix /= 0)read(line(1:ix),*)element1
     call string_trim(line(ix+1:),line, ix)
     if(ix /= 0)read(line(1:ix),*)element2
     print *,' element1 = ', element1,' element2 = ', element2
     call qp_tune(lat, qpx,qpy,ok, element1, element2)
     if(.not. ok)print *, ' Qp_tune failed'

   elseif(line(1:3) == 'QZ_')then
    call string_trim(line(ix+1:), line, ix)
    read(line(1:ix),*)Qz
    call  set_z_tune(lat%branch(0), Qz * twopi)
    print '(a20,f10.3)',' Synchrotron tune = ', lat%z%tune/twopi
   elseif(line(1:4) == 'BEGI')then
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)then
       print *,' BEGINNING <betax> <betay> <alphax> <alphay> <etax> <etay> ', &
                ' <etapx> <etapy> '
       return
     endif
     read(line(1:ix),*)lat%ele(0)%a%beta
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%b%beta
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return 
     read(line(1:ix),*)lat%ele(0)%a%alpha
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%b%alpha
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%a%eta
     lat%ele(0)%x%eta = lat%ele(0)%a%eta
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%b%eta
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%a%etap
     lat%ele(0)%x%etap = lat%ele(0)%a%etap 
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)return
     read(line(1:ix),*)lat%ele(0)%b%etap

   elseif(line(1:3) == 'ALL')then
     call string_trim(line(ix+1:),line,ix)
     read(line(ix+1:),*)value
      print '(1x,a16,1x,a11,a12,a12)','  Element name  ',' Attribute ','  New Value ','  Old Value '
     do n=1,lat%n_ele_max
      attrib=-1
      if(line(1:3) == 'SEX' .and. lat%ele(n)%name(1:3) == 'SEX')attrib=k2$
      if(line(1:3) == 'V_S' .and. lat%ele(n)%name(1:5) == 'V_SEP')attrib=vkick$
      if(line(1:3) == 'H_S' .and. lat%ele(n)%name(1:5) == 'H_SEP')attrib=hkick$
      if(attrib >0)then
        lat%ele(n)%value(attrib) = value * lat_0%ele(n)%value(attrib)
        print '(1x,a16,4x,i3,4x,2e12.4)', lat%ele(n)%name, attrib, lat%ele(n)%value(attrib), &
                                               lat_0%ele(n)%value(attrib)
      endif
      end do
     call lat_make_mat6(lat, -1)

   elseif(line(1:5) == 'LRBBI')then
     call string_trim(line(ix+1:), line, ix)
     if(ix == 0)then
       print *,' LRBBI <particle> <i_train> <j_car> <n_trains_tot> <n_cars> <current> '
       return
     endif

   elseif(line(1:3) == 'ENE')then
    call string_trim(line(ix+1:), line, ix)
    if(ix == 0)then
       print '(a,e12.4)',' Energy = ', lat%ele(0)%value(e_tot$)
     else    
      read(line(1:ix),*)energy
      lat%ele(0:lat%n_ele_max)%value(e_tot$) = energy
      print '(a,e12.4)',' Energy = ', lat%ele(0)%value(e_tot$)
      call lat_make_mat6(lat,-1)
   endif
   return
   

     if(lat%param%ixx /= 99)then
       read(line(1:),*)particle, i_train, j_car, n_trains_tot, n_cars, current
      else
       read(line(1:),*)current
     endif
      print '(1x,a12,i3,1x,a10,f10.3,a4)',' Particle = ',particle,' current= ',current,' ma'
      print '(1x,a12,i3,1x,a12,i3)',      ' i_train  = ',i_train,'    j_car = ',j_car
      print '(1x,a12,i3,1x,a12,i3)',      'tot_trains= ',n_trains_tot,'   n_cars = ', n_cars
      ans = '?'
      print '(a,$)',' Ok ?'
      read(5,*)ans
     if(ans(1:1) == 'n' .or. ans(1:1) == 'N')return


      if(lat%param%ixx /= 99)then
      print *,' Insert lrbbi for ',n_trains_tot,' and ', n_cars, ' at ' ,current,' mA'
      lat_in = lat
      call lrbbi_setup (lat_in, lat_out, particle, i_train, j_car, n_trains_tot, n_cars, current, rec_taylor)
      lat = lat_out
      lat%param%ixx = 99

     else
      lat%param%n_part = current*0.001 *(lat%param%total_length/c_light)/e_charge  
      print *,' lrbbi already inserted. Change current to ',current,' mA'  
     endif

    else
! first word is the lat element
!    element_name = line(1:ix)
    call lat_ele_locator(line(1:ix),lat,eles,n_loc,err)

     if(n_loc == 0 .or. err)then
      print *,' Cannot find element ',line(1:ix), ' Try again'
      return
    endif

    call string_trim(line(ix+1:), line, ix)
!second word is = or delta or proportional
    if(line(1:1) .eq. '=')change=setequal$
    if(line(1:1) .eq. 'D')change=delta$
    if(line(1:1) .eq. 'P')change=proportional$
  
    call string_trim(line(ix+1:), line, ix)
!third word is attribute
     attribute = line(1:ix)

     call pointer_to_attribute(eles(1)%ele, trim(attribute), .false., a_ptr, err_flag, .false.)    
     
!     if(ix_attr .eq. 0)then
     if(err_flag)then
       call type_ele (eles(1)%ele, .true., 0, .false.,0)
       print *,' Attributes of ',ele%name,' shown above. Try again'
       return
     endif

     call string_trim(line(ix+1:), line, ix)
!fourth word is value
     read(line(1:ix),*, err=10)value
     print '(a,es12.4)','value = ',value

     if(change == setequal$)a_ptr%r = value
     if(change == delta$)a_ptr%r = a_ptr%r + value
     if(change == proportional$)a_ptr%r= a_ptr%r * value
     val = a_ptr%r
     call lattice_bookkeeper(lat, err_flag)

! find slaves
     if(eles(1)%ele%key == overlay_lord$)then
       do i=1,eles(1)%ele%n_slave
         slave => pointer_to_slave(eles(1)%ele, i)
         call lat_ele_locator(slave%name,lat,slave_eles,n_loc,err)
         ix_attr = attribute_index(lat%ele(n_loc), trim(attribute)) 
         print '(a16,a10,a16,a11,i3,1x,a16,a3,es12.4)',slave_eles(1)%ele%name, &
            'slave to ',eles(1)%ele%name, ' ix_attr = ',ix_attr,attribute,' = ', val
       end do
      else
         ix_attr = attribute_index(eles(1)%ele, trim(attribute)) 
         if(ix_attr <= size(eles(1)%ele%value) .and. ix_attr > 0)then
           print '(a16,a11,i3,1x,a16,a3,es12.4)', eles(1)%ele%name,' ix_attr = ',ix_attr,&
            attribute,' = ',eles(1)%ele%value(ix_attr)
          else
           print '(a16,a11,i3,1x,a16,a3,a)', eles(1)%ele%name,' ix_attr = ',ix_attr,&
            attribute,' = ',trim(attribute)
         endif
    endif   

     
     if(eles(1)%ele%lord_status == not_a_lord$)then
         call update_ele(a_ptr%r, lat, eles(1)%ele%ix_ele, ix_attr, tot)  !set all elements with same name to same value

         ix_ele = eles(1)%ele%ix_ele
         if(ix_ele > 0 .and. ix_ele < lat%n_ele_max)then           
          print 1,lat%ele(ix_ele)%name, attribute, a_ptr%r,tot
1         format(1x, 'Element = ',a16, ' attribute ',a7,' changed  to ', f12.6, &
              ' and ',i5,' others')


          print '(a,es12.4)',lat%ele(ix_ele)%name, lat%ele(ix_ele)%value(vkick$)

          if(index(lat%ele(ix_ele)%name, 'RAW') /= 0)then
            do k = 1, lat%ele(ix_ele)%n_slave
               slave => pointer_to_slave(lat%ele(ix_ele), k, ctl)
               print *, slave%name, lat%ele(j)%value(ctl%ix_attrib)
            end do
          endif
         endif
     endif
        call lat_make_mat6(lat, -1)
     return
10   print *, ' Bad value '
 
  endif

     return
 
  end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine type_v15(ele)

  use bmad_struct
  use bmad_interface
  use bmadz_interface

  implicit none

  type(ele_struct) ele

 integer i

 if(associated(ele%mode3))then
  print '(a)', 'V matrix'
  do i=1,6
  print '(6es12.4)',ele%mode3%v(i,1:6)
  end do
else
 print '(a)',' V matrix not computed for overlay'
endif

 return
 end