anchor.diff:  					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...
  +					! perform linear interpolation between LUT values
   					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))
  
positive.code: 			v_values(i)=(i-1)/real(n_V_values-1) * (vmax-vmin) + vmin
  		enddo
  		
  		allocate(u_LUT(n_U_values,n_V_values,nx+1,nz,ny))
  		allocate(v_LUT(n_U_values,n_V_values,nx,nz,ny+1))
  		
  		write(*,*) "Percent Completed:"
  		do i=1,n_U_values
  			write(*,*) i/real(n_U_values)*100," %"
  			do j=1,n_V_values
  				domain%u=u_values(i)
  				domain%v=v_values(j)
  				
  				call linear_winds(domain,6e-5, 0, reverse_flag=.False.,useDensity=.False.,debug=.False.,fixedU=u_values(i),fixedV=v_values(j))
  				u_LUT(i,j,:,:,:)=domain%u-u_values(i)
  				v_LUT(i,j,:,:,:)=domain%v-v_values(i)
  			end do
  		end do
  		domain%u=savedU
  
negative.code: 		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,options)
  			endif
  		endif
          
  