anchor.diff:              end do
           endif
   
  -        ! if (present(debug)) then
  -        !     if (debug) then
  -        !         allocate(real_terrain(nx,ny))
  -        !         real_terrain=buffer_topo
  -        !         call io_write2d("buffered_terrain.nc","data",real_terrain)
  -        !         deallocate(real_terrain)
  -        !     endif
  -        ! endif
  -
       end subroutine add_buffer_topo
   
  +    !>----------------------------------------------------------
  +    !! Compute look up tables for all combinations of U, V, and Nsq
  +    !!
  +    !!----------------------------------------------------------
       subroutine initialize_spatial_winds(domain,options,reverse,useDensity)
  -        ! compute look up tables for all combinations of  N different U speeds and N different V speeds
           implicit none
           class(linearizable_type),intent(inout)::domain
           type(options_type), intent(in) :: options
  
positive.code:         !$omp end do
          !$omp barrier
          !$omp do
          do k=1,ny+1
              do j=1,nz
                  do i=1,nx+1
                      uk=min(k,ny)
                      vi=min(i,nx)
  
                      dpos=1
                      curdir=calc_direction(domain%u(i,j,uk), domain%v(vi,j,k))
                      do step=1,n_dir_values
                          if (curdir>dir_values(step)) then
                              dpos=step
                          endif
                      end do
  
                      spos=1
                      curspd=calc_speed(domain%u(i,j,uk), domain%v(vi,j,k))
                      do step=1,n_spd_values
                          if (curspd>spd_values(step)) then
                              spos=step
                          endif
                      end do
  
                      west  = max(i-winsz,1)
                      east  = min(i+winsz,nx)
                      bottom= max(j-winsz,1)
                      top   = min(j+winsz,nz)
                      south = max(k-winsz,1)
                      north = min(k+winsz,ny)
                      n = (((east-west)+1) * ((top-bottom)+1) * ((north-south)+1))
  
                      curnsq = sum(log(domain%nsquared(west:east, bottom:top,south:north)))
                      curnsq = curnsq / n
  
                      npos=1
                      do step=1,n_nsq_values
                          if (curnsq > nsq_values(step)) then
                              npos = step
                          endif
                      end do
  
                      top = min(j+1, nz)
                      if (top == j) then
                          bottom = j-1
                      else
                          bottom = j
                      endif
  
                      ! calculate the weights and the "next" u/v position
                      ! "next" usually = pos+1 but for edge cases next = 1 or n
                      dweight = calc_weight(dir_values, dpos,nextd,curdir)
                      sweight = calc_weight(spd_values, spos,nexts,curspd)
                      nweight = calc_weight(nsq_values, npos,nextn,curnsq)
                      ! perform linear interpolation between LUT values
                      if (k<=ny) then
                          wind_first =      nweight  * (dweight * u_LUT(spos,dpos,npos,i,j,k)  + (1-dweight) * u_LUT(spos,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * u_LUT(spos,dpos,nextn,i,j,k) + (1-dweight) * u_LUT(spos,nextd,nextn,i,j,k))
  
                          wind_second=      nweight  * (dweight * u_LUT(nexts,dpos,npos,i,j,k) + (1-dweight) * u_LUT(nexts,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * u_LUT(nexts,dpos,nextn,i,j,k)+ (1-dweight) * u_LUT(nexts,nextd,nextn,i,j,k))
  
                          u_perturbation(i,j,k)=u_perturbation(i,j,k) * (1-linear_update_fraction) &
                                      + linear_update_fraction * (sweight*wind_first + (1-sweight)*wind_second)
  
                          domain%u(i,j,k) = domain%u(i,j,k) + u_perturbation(i,j,k)
                      endif
                      if (i<=nx) then
                          wind_first =      nweight  * (dweight * v_LUT(spos,dpos,npos,i,j,k)  + (1-dweight) * v_LUT(spos,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * v_LUT(spos,dpos,nextn,i,j,k) + (1-dweight) * v_LUT(spos,nextd,nextn,i,j,k))
  
                          wind_second=      nweight  * (dweight * v_LUT(nexts,dpos,npos,i,j,k) + (1-dweight) * v_LUT(nexts,nextd,npos,i,j,k))    &
                                      +  (1-nweight) * (dweight * v_LUT(nexts,dpos,nextn,i,j,k)+ (1-dweight) * v_LUT(nexts,nextd,nextn,i,j,k))
  
                          v_perturbation(i,j,k)=v_perturbation(i,j,k) * (1-linear_update_fraction) &
                                      + linear_update_fraction * (sweight*wind_first + (1-sweight)*wind_second)
  
                          domain%v(i,j,k) = domain%v(i,j,k) + v_perturbation(i,j,k)
  
negative.code:     end subroutine initialize_spatial_winds
  
  
      !>----------------------------------------------------------
      !! Calculate the weights between the positions bestpos and nextpos
      !! based on the distance between match and indata(nextpos) (normalized by nextpos - bestpos)
      !! assumes indata is monotonically increasing,
      !! bestpos must be set prior to entry
      !! nextpos is calculated internally (either 1, bestpos+1, or n)
      !!
      !!----------------------------------------------------------
      function calc_weight(indata, bestpos, nextpos, match) result(weight)
          implicit none
          real :: weight
          real, dimension(:), intent(in) :: indata
  