anchor.diff:  		if (linear_contribution/=1) then
   			write(*,*) "Using a fraction of the linear perturbation:",linear_contribution
   		endif
  +
  +		! set up linear_mask variable
  +		nx=size(domain%terrain,1)
  +		ny=size(domain%terrain,2)
  +		if ( (.not.reverse) .and. (.not.allocated(linear_mask)) ) then
  +			allocate(linear_mask(nx,ny))
  +			linear_mask=linear_contribution
  +	
  +			if (options%linear_mask) then
  +				print*, "Reading Linear Mask"
  +				print*, "  from file: "//trim(options%linear_mask_file)
  +				print*, "  varname: "//trim(options%linear_mask_var)
  +				call io_read2d(options%linear_mask_file,options%linear_mask_var,domain%linear_mask)
  +				linear_mask=(1-(1-domain%linear_mask)*0.8) * linear_contribution
  +			endif
  +		endif
   		
   		if (options%spatial_linear_fields) then
   			if ((.not.allocated(hi_u_LUT) .and. (.not.reverse)) .or. ((.not.allocated(rev_u_LUT)) .and. reverse)) then
   				
  -				nx=size(domain%terrain,1)+1
  -				ny=size(domain%terrain,2)+1
  -				allocate(linear_mask(nx,ny))
  -				linear_mask=1.0
  -				if (.not.reverse) then
  -					if (options%linear_mask) then
  -						print*, "Reading Linear Mask"
  -						print*, "  from file: "//trim(options%linear_mask_file)
  -						print*, "  varname: "//trim(options%linear_mask_var)
  -						call io_read2d(options%linear_mask_file,options%linear_mask_var,domain%linear_mask)
  -						linear_mask(:nx-1,:ny-1)=1-(1-domain%linear_mask)*0.8
  -						linear_mask(nx,:ny-1)=linear_mask(nx-1,:ny-1)
  -						linear_mask(:nx-1,ny)=linear_mask(:nx-1,ny-1)
  -						linear_mask(nx,ny)=linear_mask(nx-1,ny-1)
  -					endif
  -				endif
  -			
   				write(*,*) "Generating a spatially variable linear perturbation look up table"
   				call initialize_spatial_winds(domain,options,reverse,useDensity)
  -				deallocate(linear_mask)
   			else
   				write(*,*) "Skipping spatial wind field for presumed domain repeat"
   			endif
  
positive.code: 		allocate(V_layers(nz))
  		allocate(preU_layers(nz))
  		allocate(preV_layers(nz))
  		
  		realnx=size(domain%z,1)
  		realnx_u=size(domain%u,1)
  		realny=size(domain%z,3)
  		realny_v=size(domain%v,3)
  		staggered = (realnx/=realnx_u).and.(realny/=realny_v)
  		
  		if (present(fixedV)) then
  
negative.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))
  					v_LUT(i,j,:,k,:)=(domain%v(:,k,:)-v_values(j))
  				end do
  			end do
  		end do
  		domain%u=savedU
  