anchor.diff:  		if (options%spatial_linear_fields) then
   			if (.not.allocated(u_LUT)) then
   				write(*,*) "Generating a spatially variable linear perturbation look up table"
  -				call initialize_spatial_winds(domain)
  +				call initialize_spatial_winds(domain,options)
   			endif
   		endif
           
  
positive.code: 					end do
  					uweight=calc_weight(u_values, upos,nextu,domain%u(i,j,k))
  					vweight=calc_weight(v_values, vpos,nextv,domain%v(i,j,k))
  					! should do some linear interpolation between this step and the next one at some point but for now...
  					domain%u(i,j,k)=domain%u(i,j,k) &
  									+    vweight  * (uweight * u_LUT(upos,vpos,i,j,k)  + (1-uweight) * u_LUT(nextu,vpos,i,j,k)) &
  									+ (1-vweight) * (uweight * u_LUT(upos,nextv,i,j,k) + (1-uweight) * u_LUT(nextu,nextv,i,j,k))
  
negative.code: 		
  	end subroutine add_buffer_topo
  	
  	subroutine initialize_spatial_winds(domain,options)
  		! compute look up tables for all combinations of  N different U speeds and N different V speeds
  		implicit none
  		type(linearizable_type), intent(inout) :: domain
  		type(options_type), intent(in) :: options
  		real, allocatable, dimension(:,:,:) :: savedU, savedV
  		real :: u,v
  		integer :: nx,ny,nz,i,j
  		
  		! the domain to work over
  		nx=size(domain%u,1)-1
  		nz=size(domain%u,2)
  		ny=size(domain%u,3)
  		
  		! save the old U and V values so we can restore them
  		allocate(savedU(nx+1,nz,ny))
  		allocate(savedV(nx,nz,ny+1))
  		savedU=domain%u
  		savedV=domain%v
  		
  		! create the array of U and V values to create LUTs for
  		allocate(u_values(n_U_values))
  		allocate(v_values(n_V_values))
  		do i=1,n_U_values
  			u_values(i)=(i-1)/real(n_U_values-1) * (umax-umin) + umin
  		enddo
  