anchor.diff:          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
  +
  +        ! we only need to go through the calculations if N is variable
  +        if (options%lt_options%variable_N) then
  +                ! $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
  -                ! 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
  +                ! $omp end do
  +                ! $omp end parallel
  +        else
  +                domain%nsquared = options%lt_options%N_squared
  +        endif
   
           ! smooth array has it's own parallelization, so this probably can't go in a critical section
           if (smooth_nsq) then
  
positive.code:         ! 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)
          else
              ! Nsq = squared Brunt Vaisalla frequency (1/s) typically from dry static stability
              stability = calc_domain_stability(domain)
  
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
  