anchor.diff:          integer :: img
   
           do img = 1, num_images()
  -            associate(ims => grids(i)%ims, &
  -                      ime => grids(i)%ime, &
  -                      jms => grids(i)%jms, &
  -                      jme => grids(i)%jme  &
  +            associate(ims => grids(img)%ims, &
  +                      ime => grids(img)%ime, &
  +                      jms => grids(img)%jms, &
  +                      jme => grids(img)%jme  &
                   )
               !$omp critical
  -            LUT(k,i,j, 1:ime-ims+1, z, 1:jme-jms+1)[i] = wind(ims:ime,jms:jme)
  +            if (this_image()==1) print*, "What does it mean that ime,jme are -1 here... this needs thought"
  +            LUT(k,i,j, 1:ime-ims, z, 1:jme-jms)[img] = wind(ims:ime-1,jms:jme-1,z)
               !$omp end critical
   
               end associate
  
positive.code:                 i = ik/n_spd_values + 1
                  k = mod(ik,n_spd_values) + 1
  
                  ! do k=1, n_spd_values
  
                      ! print the current status if this is being run "interactively"
                      if (options%parameters%interactive.and.(this_image()==1)) then
                          !$omp critical (print_lock)
                          write(*,"(A,f5.1,A$)") char(13), loops_completed/real(stop_pos-start_pos+1)*100," %"
                          !$omp end critical (print_lock)
                      endif
                      ! do j=1, n_nsq_values
                          u = calc_u( dir_values(i), spd_values(k) )
                          v = calc_v( dir_values(i), spd_values(k) )
  
                          ! calculate the linear wind field for the current u and v values
                          do z=1,nz
                              ! if (reverse) then
                              !     layer_height        = domain%z(1,1,z) - domain%terrain%data_2d(1,1)
                              !     layer_height_bottom = layer_height - (options%parameters%dz_levels(z) / 2)
                              !     layer_height_top    = layer_height + (options%parameters%dz_levels(z) / 2)
                              ! else
                                  layer_height = domain%z%data_3d(1,z,1) - domain%terrain%data_2d(1,1)
                                  layer_height_bottom = layer_height - (options%parameters%dz_levels(z) / 2)
                                  layer_height_top    = layer_height + (options%parameters%dz_levels(z) / 2)
                              ! endif
  
negative.code:             if (this_image()==1) write(*,*) "Stabilities:",exp(nsq_values)
          endif
  
  
          if (reverse.or.(.not.((options%lt_options%read_LUT).and.(error==0)))) then
              ! loop over combinations of U, V, and Nsq values
              loops_completed = 0
              if (this_image()==1) write(*,*) "Percent Completed:"
              ! $omp parallel default(shared) &
              ! $omp private(i,j,k,ik,ijk, z, u,v, layer_height, layer_height_bottom, layer_height_top, temporary_u, temporary_v) &
              ! $omp firstprivate(minimum_layer_size, n_dir_values, n_spd_values, n_nsq_values, nz,nx,ny,nxu,nyv,fftnx,fftny) &
              ! $omp firstprivate(start_pos, stop_pos, total_LUT_entries)
              ! $omp threadprivate(lt_data_m) declared at the top of the module
  
              ! initialization has to happen in each thread so each thread has its own copy
              ! lt_data_m is a threadprivate variable, within initialization, there are omp critical sections for fftw calls
              call initialize_linear_theory_data(lt_data_m, fftnx, fftny, domain%dx)
              allocate(temporary_u(fftnx - buffer*2, fftny - buffer*2, nz), source=0.0)
              allocate(temporary_v(fftnx - buffer*2, fftny - buffer*2, nz), source=0.0)
  
              ! $omp do
              do ijk = start_pos, stop_pos
              ! do ik=0, n_dir_values*n_spd_values*n_nsq_values-1
                  ik = ijk / n_nsq_values ! no +1 yet because this still needs to go through another div / mod iteration to compute i and k
                  j = mod(ijk,n_nsq_values) + 1
  