anchor.diff:  				! if we are removing linear winds from a low res field, subtract u_hat v_hat instead
   				! real(real()) extracts real component of complex, then converts to a real data type (may not be necessary except for IO?)
   				if (reverse) then
  -					if (staggered) then
  -						! Forcing data need to go back to the correct position in the grid (e.g. 2:nx-1 not 1:nx-2)
  -			            domain%u(:,z,:)=domain%u(:,z,:) - &
  -							real(real( u_hat(1+buffer:realnx+buffer+1, 1+buffer:realny+buffer  ) ))*linear_contribution
  -							
  -			            domain%v(:,z,:)=domain%v(:,z,:) - &
  -							real(real( v_hat(1+buffer:realnx+buffer,   1+buffer:realny+buffer+1) ))*linear_contribution
  -					else
  -						! if we are not on a staggered grid, just get winds from inside the buffered area
  -			            domain%u(:,z,:)=domain%u(:,z,:) - &
  -							real(real( u_hat(1+buffer:realnx_u+buffer, 1+buffer:realny+buffer  ) ))*linear_contribution
  -							
  -			            domain%v(:,z,:)=domain%v(:,z,:) - &
  -							real(real( v_hat(1+buffer:realnx+buffer,   1+buffer:realny_v+buffer) ))*linear_contribution
  -					endif
  +		            domain%u(:,z,:)=domain%u(:,z,:) - &
  +						real(real( u_hat(1+buffer:realnx_u+buffer, 1+buffer:realny+buffer  ) ))*linear_contribution
  +						
  +		            domain%v(:,z,:)=domain%v(:,z,:) - &
  +						real(real( v_hat(1+buffer:realnx+buffer,   1+buffer:realny_v+buffer) ))*linear_contribution
   				else
  -					! the internal domain is staggered by definition
  +					! note, linear_contribution component comes from the linear mask applied above on the mass grid
   		            domain%u(:,z,:)=domain%u(:,z,:) + &
  -						real(real(u_hat(1+buffer:realnx_u+buffer,1+buffer:realny+buffer  ) ))*linear_contribution
  +						real(real(u_hat(1+buffer:realnx_u+buffer,1+buffer:realny+buffer  ) ))
   		            domain%v(:,z,:)=domain%v(:,z,:) + &
  -						real(real(v_hat(1+buffer:realnx+buffer,  1+buffer:realny_v+buffer) ))*linear_contribution
  +						real(real(v_hat(1+buffer:realnx+buffer,  1+buffer:realny_v+buffer) ))
   				endif
   			
   				if (present(debug).and.(z==1))then
  
positive.code: 				call linear_winds(domain,N_squared, 0, reverse,useDensity,debug=debug)!,fixedU=u_values(i),fixedV=v_values(j))
  				debug=.False.
  				do k=1,nz
  					u_LUT(i,j,:,k,:)=(domain%u(:,k,:)-u_values(i)) * linear_mask(:,:ny)
  					v_LUT(i,j,:,k,:)=(domain%v(:,k,:)-v_values(j)) * linear_mask(:nx,:)
  				end do
  ! 				call io_write3d("v_lut"//trim(str(i))//"_"//trim(str(j))//".nc","data",domain%v)
  ! 				call io_write3d("u_lut"//trim(str(i))//"_"//trim(str(j))//".nc","data",domain%u)
  			end do
  		end do
  		domain%u=savedU
  
negative.code:         complex,parameter :: j= (0,1)
          real::gain,offset !used in setting up k and l arrays
          integer::nx,ny,nz,i,midpoint,z,realnx,realny,realnx_u,realny_v,x,y,bottom,top
  		logical :: staggered, zaxis_is_third
          real::U,V
  		real,dimension(:),allocatable::U_layers,V_layers,preU_layers,preV_layers
  		real,parameter::pi=3.1415927
  