anchor.diff:                  i = ik/n_spd_values + 1
                   k = mod(ik,n_spd_values) + 1
   
  +                print*, this_image(), ijk, i, j, k
                   ! 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 (options%parameters%interactive) then !.and.(this_image()==1)) then
  +                                !$omp critical (print_lock)
  +                                write(*,"(I5,f5.1,A)") this_image(), loops_completed/real(nz*(stop_pos-start_pos+1))*100," %"
  +                                !$omp end critical (print_lock)
  +                            endif
                               ! 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 = domain%z%data_3d(ims,z,jms) - domain%terrain%data_2d(ims,jms)
                                   layer_height_bottom = layer_height - (options%parameters%dz_levels(z) / 2)
                                   layer_height_top    = layer_height + (options%parameters%dz_levels(z) / 2)
                               ! endif
  
positive.code: 
          ! Allocate the (LARGE) look up tables for both U and V
          if (.not.options%lt_options%read_LUT) then
              allocate(hi_u_LUT(n_spd_values, n_dir_values, n_nsq_values, nxu, nz, ny)[*])
              allocate(hi_v_LUT(n_spd_values, n_dir_values, n_nsq_values, nx,  nz, nyv)[*])
              error=0
          else
              print*, "    Reading LUT from file: ", trim(options%lt_options%u_LUT_Filename)
              error=1
              ! error = read_LUT(options%lt_options%u_LUT_Filename//str(this_image())//".nc", hi_u_LUT, hi_v_LUT, options%parameters%dz_levels, LUT_dims, options%lt_options)
              if (error/=0) then
                  if (this_image()==1) write(*,*) "WARNING: LUT on disk does not match that specified in the namelist or does not exist."
                  if (this_image()==1) write(*,*) "    LUT will be recreated"
                  if (allocated(hi_u_LUT)) deallocate(hi_u_LUT)
                  allocate(hi_u_LUT(n_spd_values, n_dir_values, n_nsq_values, nxu, nz, ny)[*])
                  if (allocated(hi_v_LUT)) deallocate(hi_v_LUT)
                  allocate(hi_v_LUT(n_spd_values, n_dir_values, n_nsq_values, nx,  nz, nyv)[*])
              endif
          endif
  
  
negative.code:         integer, dimension(3,2) :: LUT_dims
          integer :: loops_completed ! this is just used to measure progress in the LUT creation
          integer :: total_LUT_entries, ijk, start_pos, stop_pos
          integer :: ims, jms
          real, allocatable :: temporary_u(:,:,:), temporary_v(:,:,:)
  
          type(grid_t), allocatable :: u_grids(:), v_grids(:)
  
          ims = lbound(domain%z%data_3d,1)
          jms = lbound(domain%z%data_3d,3)
  
          ! the domain to work over
          nz = size(domain%u%data_3d,  2)
  