anchor.diff:          ! if we are reverseing the effects, that means we are in the low-res domain
           ! that domain does not have a spatial LUT calculated, so it can not be performed
           if (use_spatial_linear_fields)then
  -            call spatial_winds(domain,rev, vsmooth, stability_window_size)
  +            call spatial_winds(domain,options,rev, vsmooth, stability_window_size)
           else
               ! Nsq = squared Brunt Vaisalla frequency (1/s) typically from dry static stability
               stability = calc_domain_stability(domain)
  
positive.code:         endif
  
          if (reverse) print*, "WARNING using fixed nsq for linear wind removal: 3e-6"
          ! $omp parallel firstprivate(nx,nxu,ny,nyv,nz, reverse, vsmooth, winsz, using_blocked_flow), default(none), &
          ! $omp private(i,j,k,step, uk, vi, east, west, north, south, top, bottom, u1d, v1d), &
          ! $omp private(spos, dpos, npos, nexts,nextd, nextn,n, smoothz, u, v, blocked), &
          ! $omp private(wind_first, wind_second, curspd, curdir, curnsq, sweight,dweight, nweight), &
          ! $omp shared(domain, spd_values, dir_values, nsq_values, u_LUT, v_LUT, linear_mask), &
          ! $omp shared(u_perturbation, v_perturbation, linear_update_fraction, linear_contribution, nsq_calibration), &
          ! $omp shared(min_stability, max_stability, n_dir_values, n_spd_values, n_nsq_values, smooth_nsq)
          !
          ! $omp do
          do k=1,ny
              do j=1,nz
                  do i=1,nx
  
                      ! look up vsmooth gridcells up to nz at the maximum
                      top = min(j+vsmooth, nz)
                      ! if (top-j)/=vsmooth, then look down enough layers to make the window vsmooth in size
                      bottom = max(1, j - (vsmooth - (top-j)))
  
                      if (.not.reverse) then
                          domain%nsquared(i,j,k) = calc_stability(domain%th(i,bottom,k), domain%th(i,top,k),  &
                                                                  domain%pii(i,bottom,k),domain%pii(i,top,k), &
                                                                  domain%z(i,bottom,k),  domain%z(i,top,k),   &
                                                                  domain%qv(i,bottom,k), domain%qv(i,top,k),  &
                                                                  domain%cloud(i,j,k)+domain%ice(i,j,k)       &
                                                                  +domain%qrain(i,j,k)+domain%qsnow(i,j,k))
  
                          domain%nsquared(i,j,k) = max(min_stability, min(max_stability, &
                                                  domain%nsquared(i,j,k) * nsq_calibration(i,k)))
                      else
                          ! Low-res boundary condition variables will be in a different array format.  It should be
                          ! easy enough to call calc_stability after e.g. transposing z and y dimension, but some
                          ! e.g. pii will not be set in the forcing data, so this may need a little thought.
                          domain%nsquared(i,j,k) = 3e-6
                      endif
                  end do
                  ! look up table is computed in log space
                  domain%nsquared(:,j,k) = log(domain%nsquared(:,j,k))
              end do
  
              if (smooth_nsq) then
                  do j=1,nz
                      ! compute window as above.
                      top = min(j+vsmooth,nz)
                      bottom = max(1, j - (vsmooth - (top-j)) )
  
                      do smoothz = bottom, j-1
                          domain%nsquared(:,j,k) = domain%nsquared(:,j,k) + domain%nsquared(:,smoothz,k)
                      end do
                      do smoothz = j+1, top
                          domain%nsquared(:,j,k) = domain%nsquared(:,j,k) + domain%nsquared(:,smoothz,k)
                      end do
                      domain%nsquared(:,j,k) = domain%nsquared(:,j,k)/(top-bottom+1)
                  end do
              endif
          end do
          ! $omp end do
          ! $omp end parallel
  
          ! smooth array has it's own parallelization, so this probably can't go in a critical section
          if (smooth_nsq) then
  
negative.code:     !! then bilinearly interpolate the nearest LUT values for that points linear wind field
      !!
      !!----------------------------------------------------------
      subroutine spatial_winds(domain,options, reverse, vsmooth, winsz)
          implicit none
          class(linearizable_type),intent(inout)::domain
          type(options_type), intent(in) :: options
          logical, intent(in) :: reverse
          integer, intent(in) :: vsmooth
          integer, intent(in) :: winsz
  