anchor.diff:                              ! v = sum( domain%v(west:east,j,south:north) ) / n
                           else
                               ! smooth the winds vertically first
  -                            u = domain%u(i,j,uk)
  -                            v = domain%v(vi,j,k)
  +                            u = u1d(i)
  +                            v = v1d(i)
  +                            ! u = domain%u(i,j,uk)
  +                            ! v = domain%v(vi,j,k)
                               ! WARNING for now domain%u,v are updated inside the loop, so the code below will end
                               ! up using u and v values that have had the linear theory applied already (not good)
                               ! eventually this should be pulled out and only the pertubration should be updated in the
  
positive.code: 
          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), &
          !$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
  
negative.code:             preU_layers(i)=sum(domain%u(:realnx_u,i,:))/(realnx_u * realny)
              preV_layers(i)=sum(domain%v(:,i,:realny_v))/(realnx * realny_v)
          enddo
          ! do i=1,nz
          !     bottom = max(1, i-vsmooth)
          !     top    = min(i+vsmooth, nz)
          !
          !     U_layers(i)=sum(preU_layers(bottom:top))/(top-bottom+1)
          !     V_layers(i)=sum(preV_layers(bottom:top))/(top-bottom+1)
          ! enddo
          U_layers = sum(preU_layers(1:nz))/nz
          V_layers = sum(preV_layers(1:nz))/nz
  
          if (.not.data_allocated) then
              if (debug) then
  