anchor.diff:  
       end subroutine initialize_spatial_winds
   
  -! Code to test that the direction calculation works.
  -!     if ((k==n_spd_values).and.(j==1)) then
  -!         print*, "Direction:", &
  -!             dir_values(i), "Calc=", calc_direction(domain%u(1,1,1), domain%v(1,1,1)), &
  -!             "  U=",domain%u(1,1,1), "  V=",domain%v(1,1,1)
  -!         if (abs(dir_values(i) - calc_direction(domain%u(1,1,1), domain%v(1,1,1)))>1e-5) then
  -!             print*, "ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR ERROR "
  -!             print*, "Error with direction:", &
  -!                 dir_values(i), "Calc=", calc_direction(domain%u(1,1,1), domain%v(1,1,1)), &
  -!                 "  U=",domain%u(1,1,1), "  V=",domain%v(1,1,1)
  -!         endif
  -!     endif
  -
   
       function calc_weight(indata, bestpos, nextpos, match) result(weight)
           ! simply calculate the weights between the positions bestpos and nextpos
  
positive.code:                         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))    &
  
negative.code:         original_buffer = options%lt_options%buffer
          variable_N = options%lt_options%variable_N
          smooth_nsq = options%lt_options%smooth_nsq
  
          stability_window_size = options%lt_options%stability_window_size        ! window to average nsq over
          max_stability = options%lt_options%max_stability                        ! limits on the calculated Brunt Vaisala Frequency
  