anchor.diff:                  end do
               endif
           end do
  -        !$omp end do
  -        !$omp end parallel
  +        ! $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
               call smooth_array(domain%nsquared, winsz, ydim=3)
           endif
   
  -        !$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 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)
           allocate(u1d(nxu), v1d(nxu))
  -        !$omp do
  +        ! $omp do
           do k=1, nyv
   
               uk = min(k,ny)
  
positive.code:     logical :: smooth_nsq
      logical :: using_blocked_flow
      real    :: N_squared
  
      !! unfortunately these have to be allocated every call because we could be calling on both the high-res
      !! domain and the low res domain (to "remove" the linear winds)
  
negative.code:                         sweight = calc_weight(spd_values, spos, nexts, curspd)
                          nweight = calc_weight(nsq_values, npos, nextn, curnsq)
  
                          ! if ((k==25) .and. (j==25) .and. (i>105) .and.(i<115)) then
                          if ((k==25) .and. (i==120).and.(j<10)) then
                              ! print*, "--------------------------------------------------------"
                              print*, j, npos, curnsq, domain%nsquared(i,j,k)
                              ! print*, dpos, spos, npos
                              ! print*, dweight, sweight, nweight
                              ! print*, " nsq=", domain%nsquared(i,j,k)
                              ! print*, curdir, curspd, curnsq
                          endif
  
                          ! 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))   &
  